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INTRODUCTION 

The  analysis  of  electromagnetic  fields  in  uniform  media  reduces  to  the  solution  of  the  vector 
Helmholtz  equation.  This  equation  is  a  vector  partial  differential  equation  of  the  form 
V  x  V  x  A  -  k2 A  =  0,  and  generally  leads  to  coupled  equations  for  the  field  components,  but  it  can 
be  separated  in  the  spherical  and  cylindrical  coordinate  systems.1  For  cylindrical  coordinates, 
detailed  treatments  can  be  found  for  localized  sources  in  the  presence  of  a  layered  cylinder,2  and 
for  localized  sources  outside  a  perfecdy  conducting  wedge,3  as  well  as  less  detailed  treatments  of 
the  conducting  wedge.4  5  In  the  present  work,  a  detailed  analysis  is  given  for  the  electromagnetic 
fields  inside  and  outside  a  uniform  sphere  immersed  in  a  uniform  medium,  in  the  presence  of  a 
localized  source  either  interior  or  exterior  to  the  sphere. 

The  solution  is  given  in  terms  of  a  vector  spherical  harmonic  expansion,  with  expansion 
coefficients  for  the  scattered  fields  expressed  in  terms  of  the  coefficients  appropriate  to  the  source 
distribution.  The  source  is  treated  using  a  clever  technique  due  to  Jackson.  A  separate  treatment 
of  the  case  of  a  perfecdy  conducting  sphere  also  is  included. 

For  explicit  results,  the  source  is  taken  to  be  a  magnetic  or  current  dipole.  Without  loss  of 
generality,  the  dipoles  are  placed  on  the  z-axis,  and  oriented  radially  or  tangentially  to  the  sphere, 
which  is  centered  on  the  origin.  Detailed  expressions  are  given  for  the  electric  and  magnetic  field 
components,  both  interior  and  exterior  to  the  sphere,  for  both  interior  and  exterior  dipoles.  These 
expressions  are  relatively  simple,  and  do  not  involve  any  recursion  formulas  for  the  coefficients, 
contrary  to  what  has  been  reported  elsewhere.7  For  completeness,  the  dc  limits  of  the  expressions 
also  are  given.  This  is  useful  for  direct  applications,  to  provide  more  transparent  forms  to  check 
for  physical  consistency  of  the  field  expressions,  and  to  provide  a  numerical  check  against  the  codes 
for  the  fields  in  the  low-frequency  limit. 

Numerical  codes  written  in  Hewlett-Packard  BASIC,  Version  3.0,  are  given  for  both  the  ac 
and  dc  field  expressions  for  the  single  case  of  exterior  dipole  and  field  point.  The  bulk  of  the  code 
given,  in  particular  the  geometrical  transformations  and  the  Bessel  function  and  Legendre 
polynomial  routines,  is  applicable  to  the  other  three  cases.  The  numerical  treatment  of  the  spherical 
Bessel  functions  incorporates  several  representations  found  in  Abramowitz  and  Stegun.8 

BACKGROUND 

For  a  time-harmonic  source  Maxwell’s  equations  for  a  uniform  medium  take  the  form 


V  x  E  =  icoB 

(1) 

VxB  =  p((o-icoe)E  +  VxMs  +  Js)  =  pJ 

(2) 

V  •  (eE)  =  p 

(3) 

V*  B  =  0 


(4) 
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where  a,  e,  and  |i  are  the  medium  conductivity,  permittivity,  and  permeability,  respectively, 
and  Ms  and  Js  are  source  magnetization  and  current  density.  The  sources  are  included  in  a  general 
way,  but  will  be  specialized  to  point  magnetic  and  current  dipoles  for  explicit  calculations. 

The  equation  of  continuity  is 

V*J-io)p  =  0  (5) 


that,  upon  substitution  of  Equations  (2)  and  (3),  can  be  written 


E  +  - 


=  0. 


O-IClE  j 

With  the  introduction  of  the  divergenceless  vector 

Js 


E'  =  E+- 


a-z( oe 


(6) 


(7) 


that  is  identical  to  E  in  a  source- free  region.  Equations  (1)  and  (2)  can  be  written 

,  VxJ5 

V  xE  =zcoB  + - — 

a-zcoe 


(8) 


and 


VxB  =  ji((o-zo>e)E'+ V  xMj).  (9) 

With  the  introduction  of  the  complex  wave  number 


£2  =  p.(tCiXJ-G)2E) 


Equations  (8)  and  (9)  can  be  written 


V  xE'  =  / 


\ 

/ 


and 


VxB=y^E'  +  ^VxMs. 


(10) 


(11) 


(12) 


Applying  the  vector  identity  V  x  V  x  A  =  V(V  •  A)  -  V*A  to  Equations  (11)  and  (12)  then 
results  in  the  vector  equations  with  sources  given  by 


(V2  +  *2)E'  =  -zcop]  -  +  V  x  Ms] 


^  V  x  V  x  Js 


(13) 
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and 


(V2  +  k2)B  =  -|i(V  xVxMj  +  Vx  Js).  (14) 

In  a  source  free  region,  the  fields  E'  and  B,  satisfying  the  associated  homogeneous  forms  of 
Equations  (13)  and  (14),  can  be  expanded  in  terms  of  transverse  electric  and  transverse  magnetic 
parts  in  the  form 


E'  =  Vx(r'F£)  +  aVxVx(r4V)  (15) 

and 

B  =  V  x  (r^)  +  pV  x  V  x  (r^)  (16) 

where  r  is  the  radius  vector,  provided  that  *F£  and  satisfy  the  scalar  Helmholtz  equation 

V2'P£  +  iV£  =  0  (17) 

and 

V2'¥M  +  k2'¥M  =  0.  (18) 

This  may  be  verified  by  applying  the  operator  V  x  V  x  (•  •  •)  to  Equations  (15)  and  (16)  in  source  free 
regions,  and  applying  the  appropriate  vector  identities. 

The  coefficients  a  and  P  can  be  determined  by  inserting  Equations  (15)  and  (16)  in  Equation 
(1 1)  in  a  source  free  region  and  using  the  identity 


V  x  V  x  (r'P)  =  /rr'F  +  V 


(19) 


to  give  the  result 


and  P  =  . 

co 


(20) 


Equation  (19)  is  a  nontrivial  identity,  and  is  developed  in  detail  below  Equation  (69)  and  preceding 
discussion). 


The  vectors  E'  and  B  can  be  determined  from  the  projections  r  •  E'  and  r  •  B,  as  can  be  seen 
by  constructing  the  projections  using  Equations  (15)  and  (16),  applying  Equation  (19)  and  the 
identity  r  •  (W  x  r)  =  0.  This  results  in 


c-  1(0 

r  •  E  =  — 

k2 


k2r2^M  +  r^ 


dr 


9^ 


M 


dr 


(21) 
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B  =  ~i(n 


7  2  df 

k2r2yE  +  r—  y  +r—£ 
dr{  E  dr 


This  shows  that  the  projections  r  •  E'  and  r  •  B  determine  4V  and 4/£,  respectively,  and  consequently 
determine  E'  and  B  via  Equations  (15)  and  (16). 

The  projections  r  •  E'  and  r  •  B  can  be  determined  everywhere  by  means  of  the  identities 

r  •  (V2  +  lc2)E'  =  (V2  +  k2)  (r  •  E')  (23) 


r*  (V2  +  A:2)B  =  (V2  +  ^2)(r*  B)  (24) 

together  with  Equations  (13)  and  (14).  These  together  give  the  inhomogeneous  scalar  Helmholtz 
equations 


(V2+*,)(r-E')  =  r-S£  =  p£ 


(V2  +  t2)(r-B)  =  r-S£  =  p£ 


where  p£  and  pB  are  given  by 


(  VxVxJ, 

■  -toapr  •  - - - +  VxMs 

k 


p8  =  -pr  •  (V  x  V  x  Ms  +  V  x  Js). 


ANGULAR  MOMENTUM  OPERATOR  AND 
VECTOR  SPHERICAL  HARMONICS 

The  subsequent  introduction  of  spherical  boundaries  near  a  localized  source  distribution  is 
greatly  facilitated  by  the  introduction  of  the  angular  momentum  operator 

L  =  -tT  x  V  (29) 

and  the  normalized  vector  spherical  harmonics 
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By  simple  manipulation  using  the  identity  V  x  (rH*)  = 
Equations  (15)  and  (16)  can  be  written 


-r  x  V'E  +  'EY  x  r  =  -r  x  V¥, 


E'  =  -iL'¥E  +  -VxL'¥M 
k2  M 


B=-iV¥M--V*V¥E.  (32) 

Since  and  ^  are  solutions  to  the  scalar  Helmholtz  equation,  they  can  be  expanded  in 
spherical  coordinate  eigenfunctions  of  the  form 

%.m  =f(kr)Ylm(0,$)  (33) 

where  f  is  a  spherical  Bessel  function  satisfying  the  differential  equation 


d2f  2df  (l2  1(1  + 1) 


J  l  J  i  *  4 

d?+T^+[k 


and  the  Y,  m  are  spherical  harmonics  satisfying,  in  particular, 

LX=/(/+i)y,,„. 

Now  note  the  following  series  of  substitutions  resulting  in  a  very  useful  identity: 

LW(*'-)r,.«(e,«] = -ir  *  VK(*r)yt„(e,(|))] 

=  -ir  X  [/;«cr)V,,K,Jje,<»]  =!fkr)  [-ir  x  VY,  ,„(6,  W 
=  V/(/  +  l)/(tr)X,..(e,«. 

If  the  scalar  functions  are  written  as  the  expansions 

^  =  i  i  1  -£==f,(kr)ra0,<» 

'=0"  =  -'\/(i  + 1) 


=  i  I  Mi‘1^Tgl(kr)Y,m(Q^) 

l  =  0m=-lyjl(l  +  \) 
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then  these  may  be  substituted  into  Equations  (31),  (32),  and  (36)  to  give 


E'=  I  I  aU*',n  x  [gtXlin] 

l  =  0m  =-i  I  k 


I  =  0  m  =  -i 


B=I  £ J a^mg,Xtn x [f,Xln] [ . 


These  forms  will  be  used  outside  of  the  source  region  for  both  the  primary  fields  due  to  the 
source,  and  the  scattered  fields  due  to  spherical  surfaces.  The  primary  fields  will  be  developed  first. 
Then,  in  a  subsequent  section,  the  scattered  fields  will  be  developed  by  applying  the  appropriate 
boundary  conditions,  and  utilizing  various  orthogonality  properties  of  the  vector  spherical 
harmonics. 

GREEN  FUNCTION  FOR  THE  PRIMARY  FIELDS 

For  a  localized  source  and  no  boundaries,  the  solutions  to  Equations  (25)  and  (26)  can  be 
expressed  in  terms  of  a  Green  function  by  means  of  Green’s  theorem: 

J  dV{G(r  -  r')  (V'2  +  k2)F(  rO -F(r')  (V'2  +  k2)G  (r  -  r')}  =  0.  (41) 

With  G  andF  satisfying 

(V'2+jfc2)G(r,r')  =  -S(r-r')  (42) 


(V'2  +  k2)F(r')  =  p(rO 


Equation  (41)  reduces  to 


F(r)  =  JdVG(r-r0p(O. 


As  shown  in  a  number  of  texts,  such  as  Jackson,7  the  free  field  Green  function  is  given  by 


Gfr-rV 


4tc|  r  -  r'| 


and  has  the  expansion  in  spherical  coordinate  eigenfunctions  given  by 

G(r-0  =  i*  S  Ukr'ri'Xkr*)  "£  « 


where  jt  and  hP  are  spherical  Bessel  functions  of  the  first  and  third  kind,  respectively,  and 
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r<l>  =  r,r  </>  r',r' otherwise. 


(47) 


EXPANSION  COEFFICIENTS  FOR  THE  PRIMARY  FIELDS 

The  primary  fields  outside  the  source  region  have  the  forms  given  in  Equations  (39)  and  (40). 
The  expansion  coefficients  can  be  determined  by  evaluating  r  •  E'  and  r  •  B.  First  note  that 


r  *  X/.m  =  0 


(48) 


and 


r  •  V/(r )  x  X,  m  =  Xl  m  •  r  x  V/(r )  =  0. 

This  leads  to  the  identity 

r  •  V  x  mr)Xi,  J  =  Mr)  r  •  V  x  Xln  =f(kr)  r  x  V  •  X,n 
=  if(kr  )L  •  X,.  =  J==LJ>',,  =  «(*r)V/(77IjV 
Thus  the  primary  field  expansion  coefficients  are  given  simply  by 


r  •  V  =  ~  i  W+ 1)  g,(kr)  Y 

k i=o 


m  =  l 


(49) 


(50) 


m  =  -l 


(51) 


and 


r  •  B'  =i  £  M Tviftkr)  "£'  a£X- 

(0i=o  --  • 


m  =-/ 


(52. 


When  Equations  (43)  and  (44)  are  applied  to  Equations  (25)  and  (26),  and  the  expansion 
Equation  (46)  is  used,  expansions  are  obtained  for  r  •  E'p  and  r  •  Bp.  Then,  Equation  (51)  becomes 


-7;  2  -W+l)g,(kr)  Y  a  £X.(9,« 

k  ‘= o  • 


m  =  / 


m  »  -/ 


(53) 


=  ik  £  z'  r,J 8,«  f  d3r'p£(r')j,(*/-<)/i,<,)(*:r=)r'„(e'> <(>'). 

For  r  <r'  for  any  source  point ,f(kr)  =  gt{kr)  =  jt(kr)  and  this  can  be  used  together  with  the 
orthonormality  of  the  spherical  harmonics  Yl  m  to  convert  Equation  (53)  to 


m,p  _ 

^ l,m 


ik 


ciW/(/  +  l) 


==  J  d,r'pE(r')hj,\krX.^A')- 


(54) 
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A  similiar  treatment  of  r  •  Bp  gives 

=^|=JdJr'pi,(r><l>«r-)l',;»(8',n  (55) 

For  points  r>r'  for  any  source  point,  hj1\krf)  is  replaced  by  jt(kr')  in  Equations  (54)  and  (55). 

APPLICATION  TO  A  SPHERE  WITH  AN  EXTERNAL  SOURCE 

For  points  interior  to  the  sphere  and  exterior  points  outside  the  source,  the  fields  have 
expansions  of  the  forms  given  in  Equations  (39)  and  (40).  In  particular,  these  forms  hold  for  all 
points  on  the  surface  of  the  sphere,  and  will  be  useful  in  imposing  boundary  conditions.  The  interior 
fields  will  be  designated  with  a  1  superscript,  and  the  exterior  scattered  fields  have  a  2  superscript. 
To  be  regular  at  the  origin,  the  interior  fields  have  the  form 


and 


m  =/ 


IC*>  MM 


/  =  0  m  =  -l 


B1  = 


oo  m  ~ 

X  X 


/  sOn  I 


aKmr)Xltm-~a?;y  x  {J^X,  J 


(56) 


(57) 


while  the  exterior  scattered  fields  are  expressed  in  terms  of  outgoing  waves  and  have  the  form 


Ww. + » [*™(V)X/,J  I 


and 


B2=  X  %jO"W)X,„-^"Vx[A»»(t2r)X,.J 

/  =  0#n  *— /  l  CO 


and,  just  outside  the  sphere  surface,  the  primary  fields  have  the  form 


~  m=l 

Ep  =  X  X 

/  —  0  in  —  — / 


af£m)x^ +jl*£y  x  imw.  j 


and 


BP  =  X  ”x'|a,M/j,(V)X,„-f  a"V x[/,(*,r)X,j[. 

CO  I 


(58) 


(59) 


(60) 


(61) 
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BOUNDARY  CONDITIONS 

The  boundary  conditions  are  continuity  of  the  tangential  E  and  tangential  H  fields  at  the 
surface  of  the  sphere.  If  R  is  the  radius  of  the  sphere,  then  the  boundary  conditions  can  be  expressed 
as 


RxE1  =  Rx(E2  +  E/>) 


(62) 


and 


— RxB'=— Rx(BJ  +  Bp). 
Hi  Hi 


(63) 


To  apply  the  boundary  conditions  in  Equations  (62)  and  (63)  it  is  necessary  to  develop  a  rather 
complex  identity.  First,  note  that 


r  x  V  x  =  r  x  I"  x  X,.„  +/(r)V  x  X,,. 

The  first  term  on  the  right-hand  side  of  Equation  (64)  has  the  form 

r  •  X,im^r-r/(r)Xl(1B  =  -r/(r)Xl>in 


(64) 


(65) 


that  is,  it  is  proportional  to  X/m.  Next,  it  will  be  demonstrated  that  the  second  term  on  the  right-hand 
side  of  Equation  (64)  also  is  proportional  to  X,  m. 

First,  note  that  by  a  standard  identity 

Vx(rxV)F  =  rV2F - (V •  r )VF  +  (VF  •  V)r - (r •  V)VF.  (66) 


The  second  term  on  the  right-hand  side  of  Equation  (66)  is  just  -3 VF.  For  the  third  term  on  the 
right-hand  side  of  Equation  (66) 


(VF  •  V)r  =  dfdfj  =  9.F5 ,  =  djF  =  VF 
and,  the  fourth  term  on  the  right-hand  side  of  Equation  (66)  is 

-(r  •  V)VF  =  -r.9,9/  =  -rfiftF 
=  -0/j  -  5^)9, F  =  VF  -  V(r  •  VF). 


Applying  these  results  to  Equation  (66)  produces  the  identity 

V  x  (r  x  V)F  =  rV2F 


-i 


_  9F  ] 
F  +  r-r-  . 
or 


(67) 

(68) 


(69) 


9 


NCSCTR  426-90 


This  result  can  be  applied  direcdy  to  V  x  Xl  m: 


v  x  Xi" = wt t x  = 7  x  (r  x  V)y'- 

-wdr^'-ii+r^}4 


v2yt.m=-^y,n=- 


Ki+ 1), 


V  x  X,  w  =  — —  rT,  +-=J='7Y, „ 
r2  V/(/  + 1) 


rxVxX-=^FTi)rxV^=-x-  <73> 

With  all  the  foregoing,  Equation  (64)  reduces  (for  spherical  Bessel  functions  /,(£/-))  to 

r  x  V  x  [/K*r)X,J  =  -[/}(&/•)  +  fcr/,(fcr)]Xim.  (74) 

Application  of  the  boundary  conditions  Equations  (62)  and  (63),  using  the  expansions 
Equations  (56)  through  (61),  and  the  identity  Equation  (74)  gives 

£  -f<  J 

Hi  Hi  Hi  J 

+i[ - 

co^ 


aSfU(W+Mj,(W} 


]X,„}  =  0 
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£  "£'  {[af:V,(*.*)  -Of Wo  -agMkJlW  x  X,,. 

I -On =— / 

r afflM*,* ) + k,Rj,(k,R))  aftW’W? )  +  kJth?\kJt)} 


CU(*.*)  +  W'iW)),„  ,  „ 

£|  JAW  “ 

The  vector  spherical  harmonics  satisfy  several  orthogonality  and  normalization  conditions 
which  are  useful  in  determining  the  expansion  coefficients6^ 

x;>,TXX,„  =  0  (77) 

Jrfnx;,„.-x„  =  8,,,8.>  (78) 

and 

JdQ(rxx;,..).(fxX,„)  (79) 

=  JdQ[(r  •  r)  (X;.„.  •  X,,)  -  (r  •  (r  •  X,„)] 

=  r%',Sn'm 

where  the  integration  is  over  all  solid  angles,  and  6 itJ  is  the  Kronecker  delta.  If  Equations  (75)  and 
(76)  are  scalar  multiplied  by  or  R  x  X*  m-  and  the  solid  angle  integral  performed.  Equation  (77) 
through  (79)  may  be  applied  to  give 

CiW )  -  rfXXkJt )  =  (80) 

CWW  + kfij,(ktR )]  -  yatWdhR  )+kJthf\kJ< )]  (81) 

=  yai,n  LMMO  "1  kjtjfikjl )] 

a^j,(k,R )  -  agh?%R )  =  )  (82) 


(*,*))  -  w^M'XkJt)  +  kjOifXkJi)] 


=  «"l WW  +  MjWO] 


where  x  =  (i,/^  and  y=  kflk%. 
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Using  the  Wronskian  {ji(u)>hil\u)}  =  i/u2,  introduce  the  determinants  of  coefficients 

Dej  =  h\\uf)  \j,iux)  +  uj/tq)]  -  [/ifW  +  UjAj’W]  (84) 

and 

DUJ  =  <  w  [/,<«,)+ « Itf’w  +  HjCM  (85) 

and  the  expressions 

ne,i  =  -y<(«2)  [//(“i)  +  uj,(u ,)]  +  u,(u 2)  +  u^Oq)]  (86) 

and 

Nm>,  =  -^//(Ma)  +  Y//(ui)  C/W  +  wJ/OO]-  (87) 


where  tq  =  /qfl  and  iq  =  •  Then,  the  solutions  for  the  expansion  coefficients,  valid  for  any 

localized  exterior  source  distribution  that  does  not  include  the  sphere  surface ,  are  given  by 


and 


IT 

u^Pe.i 


f 


Ne,i 

®E,I 


ivf 


(88) 

(89) 

(90) 


AW 

Ou,i 


(91) 


ELECTRIC  AND  MAGNETIC  FIELDS  FOR  A 
PERFECTLY  CONDUCTING  SPHERE 

For  a  highly  conducting  sphere,  the  numerical  results  using  the  coefficient  representations  in 
Equations  (84)  through  (91)  may  be  unreliable.  In  this  case  it  is  best  to  treat  the  sphere  as  a  perfect 
conductor.  Then,  the  coefficients  may  be  determined  by  means  of  the  single  boundary  condition 
Equation  (62) 


Rx(E  +  E')  =  0  (92) 

at  the  surface  of  the  sphere.  Applying  the  expansions  Equations  (58)  and  (60),  and  the  identity 
Equation  (74),  and  dropping  the  2  subscript/superscript  for  the  exterior  region,  Equation  (92) 
becomes 
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izi-‘^{a“mlh<‘XkR)  +  kRhfXkR))+a“-'mR)  +  kRj,{kRmX,,„  (93) 

I  =0  -m  K 

{afmh?XkR )  +  a". 'MR  )}RxX,.]=0. 

Again,  the  properties  of  the  vector  spherical  harmonics  are  used  in  Equations  (77)  through  (79), 
and  the  resulting  equations  solved  to  give 


M  _^M,l  M.P 

ai,m  ~  n  al,m 

uM,l 


(94) 


and 


E  _E,P 

^ l,m  i~v 

uE,l 


(95) 


again  valid  for  any  localized  exterior  source  distribution  not  containing  the  sphere  surface,  where 
now, 


and 


N^-ljfk^  +  kRjfkR)] 
Dul  =  h?\kR)  +  kRfi?XkR) 

ne.>=-WR) 

DEil  =  h?XkR). 


(96) 

(97) 

(98) 

(99) 


SPECIALIZATION  TO  DIPOLE  SOURCES 

Explicit  solutions  will  be  constructed  for  current  and  magnetic  dipole  sources.  For  these,  the 
current  density  and  magnetization  take  the  form 


Js(r)=pS(r-R0) 


(100) 


and 


Ms(r)  =  m5(r-R0).  (101) 

The  current  dipole  moment  p  can  be  specified  directly  for  a  dipole  source  shorted  to  the 
medium,  and  has  a  non- zero  value  at  dc.  To  represent  an  insulated  dipole  that  couples  only  through 
the  displacement  current,  p  should  be  written  p  =  /cop',  where  p'  is  the  eWtric  dipole  moment. 

For  V(rO  =  v8(r'  -  Ro),  with  v  a  constant  vector,  the  following  identities  are  valid: 
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V'xV  =  -vx  V'8(r'  -  Rfl) 

(102) 

V'  x  V'  x  V  =  -vV'28(r'  -  R0)  +  (v  •  V')V'5(r'  -  R^) 

(103) 

J* dr*F(T ,  r')V'S(r' -R0)  =  -V^r ,  R,) 

(104) 

J  dr*F(r ,  rOV'28(r'  -  R,)  =  V^(r ,  R*) 

(105) 

J  dr/3F(r ,  rOV'  <g>  V'8(r'  -  Ro)  =  V0  ®  V,F(r ,  R,). 

(106) 

If  Equations  (100)  and  (101)  are  substituted  in  Equations  (54)  and  (55)  and  Equations  (102) 
through  (106)  applied,  the  resulting  source  expansion  coefficients,  separated  according  to  source 
type,  are  (again,  for  an  exterior  source  distribution,  r  <r'  for  any  source  point  near  enough  to  the 
sphere  surface) 


e.cd 

^ l,m 


*  P)  *  Vo[/r/1>(^o)r;„(0o,  <J>0)] 


(107) 


and 


M,CD 

ai,m  =~ 


+ 1) 


+P*  V, 


Ica. 


} 


aE.m  =_^={iJ(Ro  .  m 


+m  •  vJ  rU%,  W^-{*A0,(W} 


(108) 


(109) 


_M,MD  _ 
ul,m  ~ 


kj\h. 

V/(/+l) 


(Ro  x  m)  • 


(110) 


There  are  two  cases  for  both  types  of  dipole,  radial  and  tangential.  For  the  radial  case 


nE,CD,R 

ul.m 


=  0 


(111) 
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^'<I,(W+F0S  w;,,<w}ly':'(e»' «  (112) 


M.MD.R  _  , 


(114) 


Without  loss  of  generality,  the  tangential  dipole  can  be  taken  to  have  only  a  9  component. 
Then  Equations  (107)  through  (110)  become 


^■^,,1WS5S''‘AW 


(115) 


M.cn.r 


1  d 


7PT)^^M(W1ae;r-(e^o) 


(116) 


g,Afc,r  1  d 


(117) 


m,md,t  _ 


m*23M2  ,i 


1  J 


(118) 


To  proceed  further,  it  is  useful  to  list  some  explicit  results  for  the  spherical  harmonics  and 
related  functions: 


s9)^- 

P"(  cos  9)  a  sin'"  9  ,  9-»0  ,  |m|£l 


(119) 


(120) 


/’.-(cose) = (-ir^~/>r(cose) 

(/  +  m)! 


(121) 


— P  "(cos  9)  =  -(/  +  m )  (/  -  m  + 1  )P"  "  ‘(cos  9)  +  m  cot  9 P  "  (cos  9). 
do 


(122) 
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Again  without  loss  of  generality,  the  dipole  may  be  taken  to  be  on  the  z-axis:  80  — » 0+ .  Then, 
application  of  Equations  (119)  and  (120)  to  Equations  (112)  and  (113)  shows  that  the  source 
coefficients  for  the  radial  dipole  vanish  unless  m  =  0,  while  application  of  Equations  (1 19)  through 
(122)  to  Equations  (115)  through  (118)  shows  that  the  source  coefficients  for  the  tangential  dipole 
vanish  unless  m  =±1. 

Before  waiting  down  the  final  results  for  the  source  coefficients,  it  is  useful  to  list  some 
special  limits: 


and 


Pf(  D=1 

P°mji±n 

sin  6  P  j  (cos  0)|  e  _^  =  — ~ 

sin  0P71  (cos  0)|  =  ~. 

1  v  e-+o+  2 


(123) 

(124) 

(125) 


(126) 


These  can  be  used  to  simplify  the  dipole  source  coefficients.  For  points  r  <  R0  the  explicit  forms 
are 


E,CD,R  n 
ai,o  =0 

(127) 

M,CD,K  pk 2^2 

(128) 

r.m.R  iunkjh  1  (2/  +  m  + 1) . 
a‘-°  -  R„  \  4*  h‘  W 

a, "*“•'=  0 

(129) 

(130) 

E,CD,T  _  E,CD,T  _ 
al,~  1  _  au 

-<aTr^h'^ 

(131) 

M.CD.T  j 

“ +  2*„  \ 

+ kJQ Sf’flyyi 

(132) 

.E.MD.T 

+  2R,  ' 

V  w + M^!"<  W 1 

(133) 
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and 


M.MD.T  M.MD.T 
al.-l  =<*1.1 


1 


irak^^  1 21  +  1 
4n 


hfXW,). 


(134) 


For  r  >R0,  h\X)(k?RQ)  is  replaced  by  j^k^Ro)  1°  Equations  (127)  through  (134). 

EXPRESSIONS  FOR  THE  SCATTERED  FIELDS 
FOR  AN  EXTERNAL  DIPOLE 

The  scattered  Fields  are  determined  by  combining  Equations  (127)  through  (134),  Equations 
(84)  through  (91)  and  Equations  (56)  through  (59)  with  m  suitably  restricted.  In  the  manipulation 
of  X,  o  and  Xz±l  some  useful  relationships  are 

2  cos  QP°,  -  sin2  0Pf  =  /(/  +  1)P,°  (135) 


and 

p 1 

— +  2  cos  0/>;  -  sin2  QPj  =  l (l  +  1  )Pj.  (136) 

sin  8 

These  appear  in  in  the  construction  of  the  radial  components  of  the  fields,  and  are  a  consequence 
of  Equation  (50). 

In  the  final  presentation  of  the  field  expressions,  only  the  Legendre  polynomials  P,  =  P°  and 
their  first  derivatives  will  appear  by  virtue  of  the  relations 

Pl 

t— -  =  —Pt  and  P!sin0  =  /(/  +  1)P, -cosG/*,  (137) 

sin  0 

and  the  Bessel  function  derivative  will  be  eliminated  via  the  relation 

f{u  )  +  uf,(u)  =  (l  +  l)Z(u)-ufl+l(u).  (138) 


MAGNETIC  DIPOLE 

.Radial 

(Interior) 

Pl 

ii 

o 

(139) 

£«=  0 

(140) 
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/-.M.ej  i  ,Wj 

9  47l/?oft  '  =  o  £>£,/ 

(141) 

s; = s  tt-c2/ + m + 1  msih^wp, 

4KRoKr  i=oDe  i 

(142) 

ge=-7^~  £(2/  +  l)?p-[(/  + 1  )jl(klr)-klrjt+l(klr)]hjlXk7R0)P, 

4nRoRr  i=o  DEl 

(143) 

B\  =0 

(144) 

(Exterior) 

© 

ii 

(145) 

o 

II 

(146) 

^si„ei(2(+  ^£,,  a;, 

9  4tcR0  i=0  ^£.1 

(147) 

4rc/V  / = o  Dej 

(148) 

=  _-m4Asing  £  p,  +  +  1)h«\t2r)-k1rh^k2r)}h"\lhRe)P, 

4n R(/  i=o  Dej 

(149) 

b]  =  0 

(150) 

Tangential 

(Interior) 


/(ony.si.e.h^  - 

47c Rr  1  =  1  D 


1 


-jt{kxr)hf\kJi0)P 


tmmii!  sin  6  -  (2/  +  1) 
4nRoRr  i  =  i/(/  +  l) 


{  tt-  K'  +  D*i®(W  -  W,<l’i<W]/\ 

uE,l 


+77“[(/  +  1);',(V)- V-/’/  +  i(^ir)J/?X1)(^?o)[/(/  +  \)P,-cos  QP ,]} 
Um.i 


(151) 

(152) 
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9  AnRoRr  ;  =  i/(/  +  l)  DMl 


(153) 


+7T-^(^'-)[(/  +  1)  Af’W-WatWl  [/(/  +  l^-cosOP,]} 


b!  =  *S1"^- £  (2/  +  l)4-[(/  +  lyi^W  -  W"’i  WV,(k,r)P, 

47 iRqRt  (  =  i  ‘-'Ej 


(154) 


(155) 


+■=■—[(/  +  lP/J-W*!^)]  [(/  +  1  [/(/  +  DP;  -COS0P,]} 

DE,l 


b; = -xinr  z  t7r->-«0^/((vW1)(W  [W + 1)/’,  -cose/-,] 

9  47c/?o/?r /  =  i /(/ +  1)  Dw/ 

+7^[(/  +  lVi(V)-Wi+i(V)]W^ 

*-'E,l 

(Exterior) 


(156) 


^  2  comii^  sin  0  sin  <|)  - 

£  = - -  2, 


I  (2/  +  W 

47tr  i  =  i 


(157) 


=  _mmM,^sm»  -  i^{_^Lir/l<'>(^)((,  +  l)h«\kjt„)-kjt^\ (k7R,)]P, 

AnRf/  /  =  i/(/  +  l)  DEi 


(158) 


+  1)P,  -COS0P,]} 

Du,i 


omfe^cosj  £  ^±^{-^K0A',,Wl('  +  iW’fVJ-^C.Wy)^/ 

9  4TC/?„r  /  =  i/(/  +  1)  £>«,/ 


(159) 


+^'-'‘,(,,(V)[('  +  '>'‘<<1,(W-W”.(Mo)]['('  +  ')P|-cos9/',]} 

UE,l 


,  /mLu/^sinOcosd)  -  A/£/  m  m 

>2 _ _ 1  v  j.  i \—E±  uV>n,  r n  i  m.OV 


47C/?oT  i  =  i 


i  (2/ + n^rc^nc/ + mr(w-WiJi(Wv,i 


(160) 


19 


NCSC  TR  426-90 


imH^cos*  -  (21  + 1)  Nkj  .  ,,,  in 

B»~  4KRor  &l(I  +  l){DM/RMik*r)h‘iWP‘ 

(161) 

+^[<;  +  mPdv)- vCi(V)]  K< + 1)CW  -  W“.  (W)  [/{/ + 1)P,  -cosep,]} 

D E,l 

=  {  ^r^c.fer)A0)(W  [/(,  +  !)/>,-  cos  6/>,] 

(162) 

+^[(i  +  i)A,(,)(/y) -VCM  [(( +  1»,0>(W  -  W“,(W]P,} 

Ue,1 

CURRENT  DIPOLE 

Radial 

(Interior) 

£,'  =  -r~t-  i  F-<2'  +  »W + 1  WV^W, 

2nklRJtri=oDu,i 

063) 

El g,m^'s“'g  J (2;  + 

2nkZRoRr  i= o  DMt/ 

064) 

II 

o 

065) 

£;=o 

(166) 

5e‘=0 

067) 

B* = -P2‘gg^  £ 1 » + 

2k RqR  <=o 

068) 

(Exterior) 

El  =  £  ^(2/  + 1)/<;  +  l)/i,mfer)A,m(y?0)P, 

AKk-Jior  i=oDMii 

(169) 

Ei=*:rr  £(^+ 

47tA:2/?or  '=o  Dm,i 

(170) 

El  =  0 

(171) 

20 


NCSCTR  426-90 


B2r  =  0 


Bl  =  0 


(172) 

(173) 


2_i>^sine  -/of  ,  niViL 


5;=  „  _  - 

■  i=o 


Tangential 

(Interior) 


I(2/  +  l)-^/I;i^rWi,(W/>, 

0  '=o 


(174) 


lottos*  £{2(  + 1)  1  m  m  + 1  (W- (WIP, 

2nkiRoRr  i= i  Aw 


(175) 


£  £l±I} {J-rRltfj,(klr)hl'XkJtlyP, 
4nkiRoRr  /«i/(/  +  l)  D£i/ 


(176) 


Aw 


[(/  +  1}/, (*/)-*, <7, „(*,r)] [(/  + 1) A,w(W- w®  (W) IW  +  1)P, -COS0P,]} 


£,‘  =  -T^  *  £  ^7Tj{^-'-«^fe<-)A,<ll(W  [((/  +  DP,  -  cosflP,] 

AidqRoRr  /=i  /(/-+- 1)  A,/ 

—  [(/+  !);,(*,  r)-V7m(*1r)][</+  mJ'to-wS’lWJP,] 

uM,l 

47t/tr  /=i  A./ 

B«  =iSr  z  §7T7TlirLry'(*1'')[('+ w,W. 

47C/CoAcr /  =  i /(/ + 1)  A,/ 

+-^-[(1  + 1)^^/)  -^i^+1(^iOW,(1)(W  [/(/  +  D^-cosOP,]} 
A,/ 


(177) 


(178) 


(179) 


7^-ry,  (*,r)[(/  +  ltfte  -  [/</  +  1«-  cose/1,]} 

uM,l 


(180) 
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(Exterior) 


,2  _  wp|i,sinBcos$  -  [  .  , ,t(1)i 


e: = -  -  — - 1  <2/ + 1  y-^-hrov) [(/ + 1 kxw  - 

AkIcJ^qT  /  =  i  £>Wi, 


(181) 


top^costj)  -  (2/  +  l)riV£,/_n 


£«  4nW  /=i/(/  +  l)lD£>< 


(182) 


+^[(/  +  l)^fer)  -  [(/  + 1  )hf\kjl0)- WffiCWl  [/(*  +  DP/  ■ -cose/*,]} 

Vm.i 


e,  =  ttyM,an»  £  (2l+_l) { ^ir^<'»(v)*m(M?o) [;(;  +  i)P,- cos6F,] 


'*  4KkJt</  ,.iHl  +  l)lDE., 


+^[((  +  l)/i,(,,(M)  -  *y*,“,(V)] «'  +  DCCMo)- MoC.WlM 

Vm,1 


(183) 


2_ip^sin0sin^  -  ,  ,^£.i,,(i) 


= — i  £  (2/ + i)^/,rw;i)(Wi 

4nr  i=i  De< 


(184) 


=  ipfe^suKD  -  gl±12 [(( +  1)Am(W _ 

4jiKor  /=i/(/  +  l)  DM>/ 

+^R  A"’(W  w  +  iv>r  W)  -  MC.(M)1  ['('  +  !)/“,-  COS  8/>,]} 
UE,I 


(185) 


jpH2*2COS<|>  (2/ +  1)  r  Ne.I  D  D  \  fff  J_  1  M.0V 


„  D  -  {“irww  to' + 1 W  W)  -  WW 

471/?,/  i=u(i  + 1)  DEii 

+7~-rhj’\k2r)  [(/  +  l)*‘,)(*2ffo) "  W®,  (Mo)l  TO  +  M,  -cos  9?,]} 

UM,l 


(186) 


DC  LIMIT  OF  THE  SCATTERED  FIELDS 
FOR  AN  EXTERNAL  DIPOLE 

The  dc  limits  of  the  field  expressions  given  in  Equations  (139)  through  (186)  are  useful  both 
to  provide  a  more  transparent  form  of  expression  for  a  consistency  check  of  the  expressions,  and 
to  provide  explicit  expressions  for  dc  applications.  These  limiting  forms  are  given  below.  Note 
that,  in  the  dc  limit,  y  -*  t5  where  5  =  CTj/CTj. 
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DC  MAGNETIC  DIPOLE 

Radial 

(Interior) 

4  =0 
4  =  0 
4  =  0 

!  Wl  y  (2l+l)l(l+l)f 

4kRqI=i  (x+l)/  +  l  |,/?oJ  1 

m^sinO  -  (2/  +  l)(/  +  l)f 

4rcR03  i-i  (X+l)/  +  l  [Ro J  ' 

4=0 

(Exterior) 

4  =  0 
4  =  0 
4=0 

m^x-l)  -  /(/  +  1)2  r £  ]+Y 

'  47t£^£  /=i(x+l)/  +  l^r  J  ' 

m^x-OsinO  -  /(/  +  !)  f RV'f R_Ylp 

8  4itf?or/?  i=i(x+l)/+l^r  J  y£0 J  ' 

4=0 

Tangential 

(Interior) 

4=0 
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(187) 

(188) 

(189) 

(190) 

(191) 

(192) 

(193) 

(194) 

(195) 

(196) 

(197) 
098) 

099) 
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Ele=  0 


£^  =  0 

mjijSinOcosiJ)  -  (2/+1)/  f  r_Y~l 

B,'~  4-rRq  /=i(x+l)/  +  lUoJ  ' 


mp.i  cos<})  -  (2/  + 1) 

e“  4tiRq  i  =  i  (x+ 1)/  + 1 


f  .  V-1 


\R°  J 


nl  m^,sin<J)-  (2/  + 1) 


4ni?o  i-i(t+l)/+l 


'  r  Y-1.. 

V^oJ  ^ 


[/(/+i)/><-cose/>,] 


(Exterior) 
£r2  =  0 

£e=0 


£2  =  0 

m^(T-l)sin9cos<l)  -  /(/  +  !)  f *Y+1(  R_Y+1£ 

"  4^/?  i=i (t+  1)/+ 1 1  r  J  J  1 


m^Ct-  l)cos<}>  ^  f RY+Y *_T+1 

6"  4tcR^  i-i(t+  1)/  + 1^  r  J  [RoJ 


[/(/  +  !)£, -cos  OP,] 


2  nrtfXzCx  —  1 )  sin  -  /  j  £ 

4jcK,>r£  1=1  (x+ 1)/  +  ll^r 


DC  CURRENT  DIPOLE 


Radial 


(Interior) 


i_  PI  "(2/ +  !)/(/  +  !) 

2710^0 1 = i  (5+l)/  +  l 


,,  _ptsin6  -  JiT"1/, 

8  2no2Ro/“1  (8+l)/  +  l  W  * 


(200) 

(201) 

(202) 

(203) 

(204) 

(205) 

(206) 

(207) 

(208) 

(209) 

(210) 


(211) 

(212) 
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El=0 
Blr  =  0 

b'6=  0 

WSrsinj  -  (2/-H)  f  r_ 

♦  2ti/?o  i-i(5+ iy+iw  ' 

(Exterior) 

p(s-D  s  /(/-n)2  r*Y+Y *  ]+1p 
'  47to2K0r/?z  =  i(5+l)/  +  ll<r  J  ^/?0  J  ' 

2_p(S-l)sine  -  fff  +  1)  f *Y+If *  Y+1* 

0  InOiRorR  z  =  i(5+l)/  +  l\,r  J  l^/?0 J  ' 

£2  =  0 

B2  =  0 

Be2  =  0 

p(S- l^sinO  -  Y+1p 

4tc£o£  ;  =  i(8+l)/  +  lVr  J  (^0  J  ' 


Tangential 

(Interior) 

i  ;usin8cos<f>  r,  (2/ -4-1)/  *6 

r~  2koJH  i-i(6+ lV  +  l^J  ' 


,  pxcos<j>  ~  (2/  +  1)  '_rY  * 

0  <=i  (5+ 1)1  + 1 


[/(/  +  !)£,-  cos  OPJ 


p x  sin <t>  “  21  +  1  f  r  Y"1 . 

♦  27tc2K03'  =  *(S+l)f  +  ll*oJ  ' 


(213) 

(214) 

(215) 

(216) 

(217) 

(218) 

(219) 

(220) 
(221) 

(222) 

(223) 

(224) 

(225) 
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„i_PHi sinO sin <})  -  2/  -4-1  (  _r_ 

r  4kRq  £i(x+\)1  +  i{r0)  1 


(226) 


2/  +  1  1 

/  +  1 

/(/  +  Dl 

.(t+l)/  +  l 

l,i,+l)P'-cosap']-^rrAi 


(227) 


—  1  _  Pj^i.  CQSfo  y  21+1  l  l  +  l 
♦  4kRq  <  =  i/(/  +  1)1(t+1)/ 

(Exterior) 


>s<}>  —  2/  + 1  [  /  +  1  .  Ixhrl  r,„  A  J  f  r  I 

l  £il(l  +  l)[(%+l)l  +  lP‘  (5+l)/  +  l[/(/+  }  '  COS0MU  (228 


2_  p(5- 1)  sinO  cos  <j>  ~  /(/  +  !)  f  £  Y+Y*  Y+Y 

'  47to2/?0Rr  i  =  i(8+l)/  +  ll/?oJ  l  r  )  1 


(229) 


2 _ p(5  1) cos<t>  “  /  f  !L 1  (H)  [/(/-}- l)/>,- cos 0/»f] 

0  47to2/?0Kr  /  =  i(8+ 1)/  + 1  {R0 J  [r  J  1 

2  p(8-l)sin0  r  <  f *TY «T*V 

♦  4nc^r  /=.(§+ l)/  +  ll/?oJ  l  r J  1 


(230) 

(231) 


p^x-VsindsinQ  -  /  +  !  f Y+Y *  Y+V 

47tRr  ,=i(x+l)/  +  ll  tf0J  l  r  J  ' 


Plijsim))  -  [  /?0(x—  1) 


4iLRRor  /  =  i[(t+l)/  +  l 


[l(l  +  l)P,-  cosOFJ- 


r(5-l) 


(6+l)/  +  lP/iUo  i  U 


/?  Y+Y*Y+1 


(232) 

(233) 


Mzii  p  _K5Z1i)_  -cose/» ilf — T+Y -T+1 

♦  4jcR/?or  (x+l)/  +  r  1  (6+l)/  +  lL  (/  )z  0P,JiUoJ  [r  J 


(234) 


APPLICATION  TO  A  SPHERE  WITH  AN  INTERNAL  SOURCE 

The  scattered  fields  still  have  the  forms  given  in  Equations  (56)  through  (59).  The  primary 
field  expansion  now  has  the  form 


E'=  £  "l  la, x[A,">(*, r)X,,J 
/=0m=-ll  k{ 


(235) 


i  T  { a**h{'\k,r )X,„  -  ~a"v  X  [f,,m(*tr)X„ 
:  =  0m  =-<  i  to 


(236) 
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The  boundary  conditions  now  are  given  by 

Rx(E1  +  E/’)  =  RxE2 
and 

— Rx(B'  +  Bp)=- RxB2. 
Hi  Ha 

In  expanded  form  Equations  (237)  and  (238)  are 


•<»  m  =  l 

I  X{| 

l -0m  =  -i 


a^h{kxR)  |  a“mph?XkxR)  a^h\\W) 

Hi  Hi  H2 


RxX, 


r  aft  WiR )  +  Wi(*i* )}  af£{h?\kxR )  +  kxRh«\kxR )} 

+i[— - — + — - 

£0Hi  WHi 


aft{/1/1>(M)  +  M^i1)(M)} 


]x/>m}=o 


and 


£  "£'  {feS,W>  +  a" */'>(*,*)  x  X,., 

/  =  0  m  =  -/ 


,  <*’ l/,(*,R ) + W^R )}  a, "/{A, ) + *,«*?>(*,«)} 

- 7; - + - *F - 

aftW’CW+tjRM’WO}™ 

- -Z  “  U' 

ki 

These  may  be  reduced  using  Equations  (77)  through  (79)  to  give 

™£V,(W ~  aftVW?)  =  -*C-  V \kxR ) 

yaft1  WiR )  +  W/(*1* )]  -  aft  V  W?  )  +  k2Rhf\k7R )] 

= 

af:lh{kxR)-af*h?XkJl)  =  -a?;'h?XkxR) 


(237) 

(238) 

(239) 


(240) 


(241) 

(242) 

(243) 
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and 


K-WW ) + *,*/,(*,*)]  -  af;M'\kJ<)+kfiil"\kJt)}  (244) 

= »(*,«) +k1Rhf\klR)] 

where  now  %  =  and  y  =  k%!kl  are  the  reciprocals  of  x  and  y.  Proceeding  as  before,  introduce 

De.i  =  [(/  +  l)ji(ui)  ~  ujl+ ,(«,)]  -;,(«,)  [(/  +  D/tfW  -  (245) 

and 

o*,,  =  Wn)  [(/  +  l)y,(«,)  -  „(«,)]  -  T/,(«|)  [(/  +  !)C  W  -  (246) 

and  the  expressions 

A^,  =  A.(‘W  [(/  + 1  )h?\uf)  -  -  tft?  W  [(/  +  1)*,%,)  -  MjA&GO)  (247) 

and 

A?*.,  -*&,%,)  [(/  + 1  )h?\uj  -  ihfcffiWl  -  W  W  t« +  ifcPW  -  (248) 


where  again,  ux  =  A:^  and  112  =  k?R.  Then,  the  solutions  for  the  expansion  coefficients,  valid  for  any 
localized  interior  source  distribution  that  does  not  include  the  sphere  surface,  are  given  by 


Ne,i 

&E,l 


(249) 


and 


ix 

u\De,i 


Ay, 1 

Dm.i 


M, 

ll,m 


(250) 

(251) 


iyf  _m,p 
«1  DMtl 


(252) 


For  dipole  sources,  the  expressions  in  Equations  (127)  through  (134)  may  be  used,  with 
hj'XkjRo)  replaced  by  j,(kxR0)  and  then  the  2  subscripts  replaced  by  1.  All  the  angular  information 
contained  in  Xi0  and  Xu,  remains  the  same. 
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The  end  result  is  that  the  field  expressions  for  the  source  interior  to  the  sphere  can  be  obtained 
directly  from  Equations  (139)  through  (186)  for  the  exterior  source  case,  by  means  of  simple 
exchanges  and  substitutions. 

To  obtain  the  new  field  components  for  region  1,  replace  ( )  by  j,  ( )  everywhere  in  the 

old  expressions  for  region  2,  replace  y,  N,  and  D  by  y,  N,  and  D  in  these  expressions,  and  then 
replace  2  by  1  everywhere  it  appears  explicitly. 

To  obtain  the  new  field  components  for  region  2,  first  define  [/,(«)]  =  (/  + 1 )f(u)  -ufui(u), 
then  replace  ji{kxr)  and  [/,(&,/-)]  by  hW^r)  and  [hWk^r)],  h\x\hji0)  and  [h^Xk^Ro)]  by  j^Ro)  and 
[/;(£i^o)]>  all  unbarred  quantities  by  barred  quantities,  change  the  remaining  subscripts  from  1  to  2, 
and  reverse  the  sign  of  DMl.  The  results  of  these  operations  are  listed  in  the  following  section. 


EXPRESSIONS  FOR  THE  SCATTERED  FIELDS 
FOR  AN  INTERNAL  DIPOLE 


MAGNETIC  DIPOLE 

Radial 

(Interior) 

^  ** 

n 

o 

(253) 

<& 

ii 

o 

(254) 

.  (omLi.L  sin  0  ~  NE . 

Ei = — tts —  r  (2/ + 1 

47l/\o  1  =  0  L)  E'i 

(255) 

.  imu.L  -  NE , 

B  S  <2' +  1>,y  + 1 

ATIKqT  i=o  De  1 

(256) 

.  /mii.L  sin0  -  NEI 

B  —  ith  Z(2l  +  l)y^LKl+im,r)-k,rj:„(klr)y,(k<Ra)P, 

47C«0r  i=o  oEi 

(257) 

Bx=  0 

(258) 

(Exterior) 

E]  =  0 

(259) 

Ee2  =  0 

(260) 
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El  =  K'7'Tr9'  * (2' +  1)=-*,“W)A( W, 

?  4^/?o/?  / =0  D£i/ 


(261) 


B*  = 


2  J-(2/ +  !)/(/+  l)*,">(fcr)y,(t1fi0)P, 


'  4nRoRr  i=oDe  t 


(262) 


mi^sinG  - 


1 


=  1  (21  + 1)=— £(/  +  D^W)" 

47C/?o«^  '=0  Z)£i, 


fl;=o 


(263) 

(264) 


Tangential 

(Interior) 

= _«M**L*  £ (2(  +  i&UjtoWWt. 

4nr  i=i  DUJ 


(265) 


_ _omyi|smj  -  gUl} {_Nej  [(( +  _ Wj  ,(*lRo)J/>, 

4jri?(/-  <«»/(/ +  1)  £>£,/ 

+^^[(/  + 1  )yi(^x'')  -  Krii  ^i(kir)]RJi(kiR0)  [1(1  +  l)P,-cosQP,]} 
Dm, i 


(266) 


♦  47Cfior  i=i/(/  +  l)  DWi/ 

+=^2-rj,(klr)  [(/  +  1)./,(M0)- W/+i(Wl  W  +  1VV -cos0P,]} 
De.i 


(267) 


in^Asmecosj^  (;  +  1  w,  ,(W1/>, 

'  4rctf0>'  t-i  £>£,i 


(268) 


=  £  Ql±l} { 

4nRor  i  =  i  1(1  +  1)  DMti 


(269) 


+^%  +  lWftr)- VAmOWI  [(/  +  mw- W,.«W  + W  -coseB,]} 

De.i 
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_  _  mi^sin$  £  ^±12  {  rRJcfj,(ktr)j,(k,R0)  [1(1  +  1  )/>,  -  cos  6>J 
f  4ji/?or  /=i/(/  +  l) 


(270) 


De.i 


(Exterior) 


_  _-'~esioj  j  +  ,,1^^ 
47t/?r  i=i  DMj 


(271) 


l2  iwmMaSin^  -  (2Z  +  1).  1  .m 


E  —  — 
C0 


.  ljfTZ^{=-rh ^rmi  +  DUk.Ro)-, WmW. 

4jt/?o/?r  /=i/(/  +  l)  Z)£i 


(272) 


+=?-[(/  +  l)C(^)-VC(W,W[/(/  +  DP, -cos  OP,]} 
Dm, i 


jornfecos'1’  £  H±l]{  1  RanklRjVf  +  r)h?Xiv)-lvh«! ,(**•)]*, 

♦  47t/?ol?r  /  =  i  /(/  +  1) 


(273) 


+ J-r^(1)(V)[(/  +  D-zXW- W/+,(W1  [/(/  +  D^-cosOP,]} 
De.i 


mfesmecos'f  £  +  1  [((  + 1  - kRj l(*lR0)]fc<‘,>(*2r JJ*, 

47cPoPr  /=i  D£-, 


(274) 


B‘ =  ^  m  +  n { =~rR’£h‘)(lvm'R°)Pl 

AiiRoRr  1=1 1(1  +  1)  £)M/ 


(275) 


-=}-[(/  +  l)k,,%r)-k2r/l"\(k2r)\  [(I +  1)a(W"  [/(/  + 1)/>,  -cos6^,]} 


=  {  =i-  rR^WdylMW  [1(1  + 1  )P,  -  cos  9/-J 

?  4^/?oPr  i=i  1(1+  1)  £>Mi 

D  r  I 


(276) 
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CURRENT  DIPOLE 

Radial 

(Interior) 


= --rrF-  £  =!=!  (2; + mi + 1  tMKDMWP, 

i=qD M  i 


{111) 


El  £  (2/  + 1)=%/  + 

47C*iR0''  '=<>  Omi 


(278) 


£,x=0 

fi!  =  0 


5e=0 


(279) 

(280) 
(281) 


T  4jik0  / =o  DUI 


(282) 


(Exterior) 


eS.'Tt--  £  =J— cm + dki + 1)  h?'ovmwi 


r 


InkfRoRr  i=oDMi 


(283) 


El  =  -impl^S1^  £  (2/  + 1)  =J-[(/  + 1  tf'(V)  -  VCWUKWIf. 


hdc^RoRr  i=o 


Dm, i 


(284) 


£j  =  0 


£‘  =  0 


fie=0 


(285) 

(286) 
(287) 


~  'Ti!  9  £(2/+i  )=-n|V)i(lAlf, 

’  2ickoK  /  =0  Dw< 


(288) 
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Tangential 

(Interior) 


(op Hi sinO cos <}>  -  Nmj 

AnkxRor  1  =  1  DU't 


j,(kxr)Ki  + 


(289) 


„i  apiLiCOsQ  -  (2/  +  1)  rNEtl  _d  , 2;/|  -w/t  D  ,6 
£9  ~  — a  b  d  "  t/i  ,  i\  itt  rR0kxjl(klr)jl(/cxK0)/Jl 
4kKiRoT  i  =  i/(/  +  1)  DEl 


[(/  +  Dy^r )  -  Vy, + 1(^)]  W  +  -  Wi + .( W1  [/(/  + 1 V*/ "  COS0PJK29O) 

D  i 


'  M,l 


=  mp|X,sinj  £  )y(W[(f(  +  i)Pj  -cose?,] 

*  4rot,F<>r  1  =  1  /(/  +  1)  DEl 


N 


+=^[( l  +  \)ji(kxr)-kxrjt+  ,(*,/■)]  [(/  +  l)j,(kxR0)  -  kxRJl+x(kxR0)]P  t} 


D 


M,l 


(291) 


_ ipftt, ^sinesin^  £  + 

4jtr  1=1  D£(, 


(292) 


jp^sinj  £  ei±i]{_^Wry(t  r)[(/  +  ly  (W- W, .,(«]/>, 

4kRoT  i  =  i /(/  +  1)  Dm  / 


(293) 


+^+RJ,(klR0)  [(/  +  lVKV) "  Wi+i(V)l  W  +  l)^i  -  cos  QP,]} 
De.i 


=  _/p^cos*  -  [(/  +  1  r/,.,(*1r)]P; 

?  47CK0r  1  =  1 /(/  +  1)  £)£/ 

+^rjl(klr)  [(/  +  lMW  -  Wi+iOW]  W  +  D^i  -  cos  OF,]} 
Dm, i 


(294) 


(Exterior) 


2  i(op^sinecos(|)  -  ,  ,,  1  ,.(i) 


Et  = 


2nkfRoR  r  i  =  i  D 


I  (2/  + 1)= — *rX^2r)[(/  +  l)j,(kxR0)  -  kxRJl+x(kxR0)]P , 


M.l 


(295) 
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2  iCPpMgCOsfr  -  (2/  +  l)f  1  -.2.(1), 


47t*12*0/?r  “,/(/  +  !)  l£>£, 


{= , 


(296) 


+ ==-[(/  +  W  +  H/iW-WuiW  [/(/  +  -cosO/*,]} 

Dm, i 

*-7&r  *  §7TT{=-'-*A2*r(V)/ff.*J['(/+ w.-cose/-,] 

*  47dkfi?oRr  i«i /(/  +  !)  D£i, 


+==-  [(/  +  DtfW)- -  W  + 1  )//0 W- -  W/ + 1(  WFJ 

D  mi 


Bt  =  /?^S1-e---  £  (2/  + 1)  J-h^l\k2r)jl(k1R0)Pl 


4k Rr  /= i  Z7£/ 


2_p^sin(J)  -  (2/  + 1)  r  2y  _,.(i) 


5e  = 


4jiRoRr  t=i  1(1 +  1)  Du,i 


{-=4-rhX\kf)  [(/  +  1)/,(W  "  WwCWi 


+ J-  [(/  + 1)  hfXkj)  -  */*&(^)Wi(W  [/(/  +  1)P,~  cos  0P,]} 
De.i 

££CQS*  -  H±i2{  ‘  RJ^IQUl  +  l)fc<‘>(V)-*2'-'>">.(^)^, 

y  4KRoKr  /  =  i  /(/  +  1)  D£i, 

-=^-  rhfXkjr)  [(/  +  1)/,(W  ■ -  WmWl  W  +  W-costf,]} 
Dm.i 


(297) 


(298) 


(299) 


(300) 


DC  LIMIT  OF  THE  SCATTERED  FIELDS 
FOR  AN  INTERNAL  DIPOLE 

The  dc  limits  of  the  field  expressions  given  in  Equations  (253)  through  (300)  are  not  simply 
derived  from  the  dc  expressions  for  the  external  dipole,  given  in  Equations  (187)  through  (234). 
The  limiting  forms  are  derived  directly  and  are  given  below.  In  writing  the  expressions,  the  results 
are  given  in  terms  of  x  =  IVH2  and  5  =  Oi/o2. 

DC  MAGNETIC  DIPOLE 

Radial 

(Interior) 
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£e=0 

e;=o 

nvtfr-l)-  l2(/  +  l)  fryfRoY 

4nRorR  ,=i(z+l)I  +  l{R  J  [r  J  ' 

.  mm(T-l)sine  -  /(/  +  !)  f  r  W ^\p 

6  4-nRj-R  ii(z+l)l  +  l{R  ){r  )  1 

B\=  0 
(Exterior) 

£2  =  0 

£e=0 
£2  =  0 

r,2 _  mfet  “  (2/  +  !)/(/ +  l)f £oY 

'  47t£oT2'  =  '  (t+lX  +  1  U  J  ' 

m^TsinO  -  (2/  +  1)/  f  . 

6  4riV2  Ad+lV  +  l^r  }  1 
£2  =  0 

Tangential 

(Interior) 

£,'  =0 
£0=0 
£^=0 

mp-tCx- l)sin9cos<j>  ~  /(/  +  !)  f  /O'* 

4nRorR  /  =  i(t+ 1)/  + 1(£  Ji,£  J  ' 


(302) 

(303) 

(304) 

(305) 

(306) 

(307) 

(308) 

(309) 

(310) 

(311) 

(312) 

(313) 

(314) 

(315) 

(316) 
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m^pt-  l)cos<})  -  (/  +  !)  ( 

B»~  4  nRorR  /=i(x+l)/  +  H 

j  l)sin<J>  -  /  +  1  ( r_ 

B*  ~  4nRorR  /  =  i  (x  + 1)/  + 1  [r 

(Exterior) 


[l(l  +  l)Pt-  cos0P,] 


E2  =  0 
Ee=0 


£*  =  0 


Y 


B; 


,  m^TsinOcoscj)  -  (2l  +  l)(l  +  l)f Bo 
4kRoT2  /=i  (T+l)/  +  l 


Pi 


n2  m^Tcost})  -  (21  + 1) 

Bq  —  *  «  i 


47t/?(/‘2  /-i(t+l)/  +  l 


v  '  y 


[/(/  +  l)P/-cosePJ 


m^xsin$  -  (21  + 1)  ( BoY  • 

♦  47iBor2  <=i(T+l)/  +  ll,r  J  ' 


DC  CURRENT  DIPOLE 
Radial 

(Interior) 

r»  E£zR  y  /2(/  +  1>  f iTf 

'  “  47ca1B(/£  (8+ 1)/  + 1 U  J  {R  )  1 


„i  p(S-  l)sin0  -  /(/  +  !) 

e  47iaji?or/?  /  - 1  (8+ 1)/  + 1 


Bo 

B 


£*=0 

Brl=0 

Be'=0 


(317) 

(318) 

(319) 

(320) 

(321) 

(322) 

(323) 

(324) 


(325) 

(326) 

(327) 

(328) 

(329) 
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i  p  (S-l)^  sin6  -  /  (  r][  *0^ 

4 kRoR  £i$+1)1  +  i{r){r  J  1 

(Exterior) 

-3  P  y  QLlMlllfM'p 

r  2it<r2fi&r2/=i  (5+l)/  +  l  (rj  ' 
psin6  -  (2/  + 1)/  ( Rtlp 

0  27tG2K(/'2'=1(8+l)/  +  lV r  J  ' 

E?  =  0 
B?  =  0 
B]  =  0 

Pl^sinO  -  (2/  +  1)  f 

*  2jcfl0r  i=i  (5+1)/  +  !^  r  )  1 


Tangential 

(Interior) 

p(5- l)sinOcos<l>  “  /(/  +  !)  ^  'p 

£r_“  4no1fl0Kr  i=i(6+l)/  +  l^  J  ' 


p(8-l)cos<l>  ~  /  + 1  f  r  Tf^°T 

8-_  4tc OiRoRr  £i@+1)1  +  i{r  )  [r  j 


[/(/  +  !)/>,- cos  OP,] 


p(5- 1 ) sin <}>  -  /  + 1  f  r_V( R^lp 

'♦  47CG^o«^  1.1  (5+ 1)/  + 1  [r  )  t  R  J  ' 


,  pp.i(T- 1)  sin  Osin  (j)  -  / 

B'~~  4nRr  i.i(x+l)/  +  l 


.i _  PM-i sin<j)  " 

9  4nRRori=i. 


*o(t-1) 

(t+1)/  +  1 


[/(/  +  l)P,-cos0/>,] 


r(S-l) 

(§+!)/  + 1  1\{r)[r) 


(330) 

(331) 

(332) 

(333) 

(334) 

(335) 

(336) 

(337) 

(338) 

(339) 

(340) 

(341) 
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m cosj)yf  /?o(x-l)  ^  r(S-l)  r/(/  +  n/>_cosep1if iTf ?2V 
5*  47iR/?oT  i?i|(T+l)/  +  lP'  (8+l)/  +  l  '  'jUJl* 


(342) 


(Exterior) 


V 


2  p  sin 8 cos (j)  “  (2/  +  !)(/  +  l)T^o 
r"  2710^  '  =  >  (8+ 1)/  + 1  (  r  , 


P, 


El  =  i  75TT77Tt[  -  W*  +  1  )P-  COS  OP  J 


27ta2/?(/2'=i(8+ 1)/  + 1 


\r ) 


_2  _  p  sin  ~  (2/  + 1) 

^  —  't  2^ 


^  2jta2/?(/'2,=1(S+l^+^ 


rM 


2  pxHjSinOsiiKj)  -  2/  +  1  f /?<> 

B,‘  4tcr2  ii(x+l)l  +  l{r 


V 


(343) 


(344) 


(345) 


(346) 


B>=£^£(2/  +  1). 


ZRr 


4nRor2  i= i 


8;=^i(2/tiv 


47c/?or  /=i 


(/  +  l)[(l+l)/  +  l] 
xR* 


[1(1  +  1)P,  -  cos  OF,]  - 


2r 


/[(8+l)/  +  l] 


>]f*v 


(347) 


V' 


2r 


l(/  +  l)[(r +!)/  +  !]  '  /[(8+l)/  +  l] 


-[/(/  +  l)P,  -  cos  QP,] 


'rA1 


(348) 
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APPENDIX  A 

NUMERICAL  METHODS  AND  HEWLETT-PACKARD  BASIC  3.0  COMPUTER 
CODES  FOR  CURRENT  AND  MAGNETIC  DIPOLE  FIELDS 

In  tills  Appendix,  codes  are  given  fcr  computing  the  electric  and  magnetic  fields  external  to 
a  sphere,  for  dipole  sources  external  to  the  sphere.  Codes  are  given  for  ac  sources(MAGNDIPSPH 
and  CURRDIPSPH)  and  for  dc  sources  (MGDIPSPHDC  and  CRDIPSPHDC).  The  Legendre 
polynomials  and  their  first  derivatives  are  generated  in  the  subroutine  Pl_pldot  based  on  the 
recursion  formulas.1*11  The  spherical  Bessel  functions  can  be  generated  from  closed-form 
expressions, A1  but  these  have  large  roundoff  errors  for  any  given  argument  as  the  order  increases. 
Instead,  the  Bessel  functions  are  generated  from  various  representations,  depending  on  the  values 
of  argument  and  order.  These  routines  were  taken  from  Abramowitz  and  Stegun,A1  and  are  discussed 
in  detail  below. 

The  structure  of  the  remainder  of  this  Appendix  is  as  follows: 

MAGNDIPSPH 

CURRDIPSPH 

MGDIPSPHDC 

CRDIPSPHDC 

The  subroutines: 

Geomdipsph 

Geomfldpos 

Geomfldbe 

Jcomb 

Hcomb 

Spherejnz 

Spherehnz 

Jnuevrywhr 

Hnuevrywhr 

Plpldot 

Ndm 

Nde 

Ndtnicon 

Ndeicon 

Discussion  of  series  representations  of  Bessel  functions. 

The  subroutines: 

Jn 

Hln 

Jnu 

Hlnu 

Prinlogz 

Atn2 

Gamma 

Discussion  of  asymptotic  representations  of  Bessel  functions  for  large  argument 
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The  subroutines: 

Jnasy 

Hlnasy 

Pnz 

Qnz 

Discussion  of  uniform  asymptotic  representation  of  Bessel  functions  for  large  order. 

The  subroutines: 

Jfiuniasym 

Hfkuniasym 

Jhlproduni 

Uk 

Ukcoefs 


PROGRAM  MAGNDIPSPH 
10  OPTION  BASE  1 

20  DIM  Xf(3),Xd(3),M(3),Rds(3,3),Bp(3,2),Ep(3,2),B(3,2)3(3,2) 

30  DIM  Pl(100)d>ld(100),Ndma(100,2)>Ndea(100^),Et(3,2)3t(3^) 

40  INPUT  "FREQUENCY  AND  MEDIUM  AND  SPHERE(Ss<0 

INF)CONDUCTTVITIES(MHO/M)?,,3,Sm,Ss 

50  INPUT  "RELATIVE  PERMEABILITY  OF  MEDIUM  AND  SPHERE?",Mm,Ms 
60  INPUT  "RELATIVE  PERMITTIVITY  OF  MEDIUM  AND  SPHERE?"3mJEs 
70  INPUT  "POSITION  OF  FIELD  POINT?(M)",Xf(*) 

80  INPUT  "POSITION  OF  SOURCE  POINT?(M),,,Xd(*) 

90  INPUT  "MAGNETIC  DIPOLE  MOMENT  VECTOR?(AMP-MA2)",M(*) 

100  INPUT  "SPHERE  RADIUS(M)  AND  MAXIMUM  POLAR  ANGLE  INDEX?",A,E11 
1 10  INPUT  "RELATIVE  ERROR  FOR  TRUNCATION?"^ 

120  Lmax=Ell+l 

130  W=2*PI*F 

140  M0=4*PI*l.E-7 

150  E0=8.85415E-12 

160  Utl=W*M0 

170  Ut2=W*E0*Utl 

180  K12r=Ut2*Es*Ms 

190  K12i=Utl*Ms*Ss 

200  K22r=Ut2*Em*Mm 

210  K22i=Utl  *Mm*Sm 

220  Md=SQR(K12r*K12r+K12i*K12i) 

230  K  lr=SQR((Md+K  1 2r)/2) 

240  Kli=SQR(ABS(Md-K12r)/2) 

250  Md=SQR(K22r*K22r+K22i*K22i) 

260  K2r=SQR((Md+K22r)/2) 

270  K2i=SQR(AB  S(Md-K22r)/2) 

280  Mefac=W*M0*Mm/4/PI  !  Electric  fields  in  volt/meter 
290  Mefacr=Mefac*K2r 
300  Mefaci=Mefac*K2i 

310  Mmfac=M0*Mm/4/PI  !  Magnetic  fields  in  tesla 
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320  Mmfacr=-Mmfac*K2i 
330  Mmfaci=Mmfac*K2r 
340  Mur=Ms/Mm 

350  !Only  frequency  and  medium  constants  to  here 
360  REDIM  Ndma(Lmax,2),Ndea(Lmax,2) 

370  IF  Ss>0  THEN  GOTO  410 

380  CALL  Ndmicon(Lmax,K2r,K2i,A,Ndma(*)) 

390  CALL  Ndeicon(Lmax,K2r,K2i,A,Ndea(*)) 

400  GOTO  430 

410  CALL  Ndro(Lmax,Klr,Kli,K2r,K2i,Mur,A>Ndma(*)) 
420  CALL  Nde(Lmax,K  lr,K  1  i,K2r,K2i,Mur,A,Ndea(*» 
430  !  Added  only  sphere  radius  to  here 

440  CALL  Geomdipsph(Xd(  *)  ,M(*),kds( * ),R0,Mr,Mt) 

450  CALL  Geomfldpos(Xf(*),Rds(*),R,Ct,St,Cp,Sp) 

460  Zr=K2r*R 

470  Zi=K2i*R 

480  Z0r=K2r*R0 

490  Z0i=K2i*R0 

500  REDIM  Pl(Lmax)  Fld(Lmax) 

510  CALL  Pl_pldot(Lmax- 1  ,CXPl(*),Pld(*)) 

520  MAT  Ep=  (0) 

530  MAT  Bp=  (0) 

540  L=2 

550  MAT  Et=  (0) 

560  MAT  Bt=  (0) 

570  CALL  Hcomb(L- 1  ,Zr£i  Jl,Hzr,Hzi,Chzr,Chzi) 

580  CALL  Hcomb(L- 1  ,ZOr,ZOi,R0,HzOr,HzOi,ChzOr,ChzOi) 
590  Te  lr=Hzr*Hz0r-Hzi*Hz0i 
600  Tel i=Hzr*Hz0i+Hzi*Hz0r 
610  Te2i^Telr*Ndma(L,l)-Teli*Ndma(L,2) 

620  Te2i=Telr*Ndma(L,2)+Teli*Ndma(L,l) 

630  Tb2r=T e  lr*Ndea(L,  1  )-Te  1  i*Ndea(L,2) 

640  Tb2i=Te  lr*Ndea(L,2)+Te  1  i*Ndea(L,  1 ) 

650  Et(  1 , 1  )=Mt*St*Sp*Te2r*Pld(L)/R 
660  Et(  1 ,2)=Mt*St*Sp*Te2i*Pld(L)/R 
670  Te3r=Hz0r*Chzr-Hz0i*Chzi 
680  Te3i=Hz0r*Chzi+Hz0i*Chzr 
690  Te4r=Te3r*Ndma(L,  l)-Te3i*Ndma(L,2) 

700  Te4i=Te3r*Ndma(L,2)+Te3i*Ndma(L,  1) 

710  Te5r=Hzr*Chz0r-Hzi*Chz0i 
720  Te5i=Hzr*Chz0i+Hzi*Chz0r 
730  Te6r=Te5r*Ndea(L,l)-Te5i*Ndea(L,2) 

740  Te6i=Te5r*Ndea(L,2)+Te5i*Ndea(L,l) 

750  Bt(  1 , 1  )=Mr*Tb2r*L*  (L- 1 )  *P1(L)/R/R0 
760  Bt(l,2)=Mr*Tb2i*L*(L-l)*Pl(L)/R/R0 
770  Bt(  1 , 1  )=Bt(  1 , 1  )+Mt*Cp*St*Te6r*Pld(L)/R 
780  Btf  1  ,2)=Bt(  1 ,2)+Mt*Cp*St*Te6i*Pld(L)/R 
790  Tb3r=Chzr*Chz0r-Chzi*Chz0i 
800  Tb3i=Chzr*Chz0i+Chzi*Chz0r 
810  Tb4r=Tb3r*Ndea(L, l)-Tb3i*Ndea(L,2) 

820  Tb4i=Tb3r*Ndea(L,2)+Tb3i*Ndea(L,  1) 

830  Tb5r=K22r*Te2r-K22i*Te2i 
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840  Tb5i=K22r*Te2i+K22i*Te2r 

850  Tb6p=Te3r*Ndea(L,  1  )-Te3i*Ndea(L,2) 

860  Tb6i=Te3r*Ndea(L,2)+Te3i*Ndea(L,l) 

870  Et(2, 1  )=Et(2, 1  )+Mt*Sp*Te4r*(Pl(L)-Ct*Pld(L)/L/(L- 1 )) 
880  Et(2,2)=Et(2,2)+Mt*Sp*Te4i*(Pl(L)-Ct*PId(L)/L/(L-l)) 
890  Erf 2, 1  )=Et(2, 1  )-Mt*Sp*Te6r*Pld(L)/L/(L- 1 ) 

900  Et(2,2)=Et(2,2)-Mt*Sp*Te6i*Pld(L)/IV(L-l) 

910  Bt(2, 1  )=Bt(2, 1  )-Mr*St*Tb6r*Pld(L)/R0 
920  Bt(2,2)=Bt(2,2)-Mr*St*Tb6i*Pld(L)/R0 
930  Bt(2, 1  )=Bt(2, 1  )+Mt*Cp*Tb5r*Pld(L)/17(L- 1 ) 

940  Bt(2,2)=Bt(2,2)+Mt*Cp*Tb5i*Pld(L)/L/CL-l) 

950  Bt(2,l)=Bt(2,l)+Mt*Cp*Tb4r*(Pl(L)-Ct*Pld(L)/L/(L-l)) 
960  Bt(2,2)=Bt(2,2)+Mt*Cp*Tb4i*(Pl(L)-Ct*Pld(L)/L/(L- 1 )) 
970  Te7r=Telr*Ndea(L,l)-Teli*Ndea(L,2) 

980  Te7i=Te  lr*Ndea(L,2)+Te  1  i*Ndea(L,  1 ) 

990  Te8r=Te3r*Ndma(L,l)-Te3i*Ndma(L,2) 

1000  Te8i=Te3r*Ndma(L,2)+Te3i*Ndma(L,  1 ) 

1010  Et(3, 1  )=Et(3, 1  )+Mt*Cp*Te6r* (Pl(L)-Ct*Pld(L)/IV(L- 1 )) 
1020  Et(3,2)=Et(3,2)+Mt*Cp*Te6i*(Pl(L)-Ct*Pld(L)/lV(L-l)) 
1030  Et(3, 1  )=Et(3, 1  )-Mt*Cp*Te8r*Pld(L)/L/(L- 1 ) 

1040  Et(3,2)=Et(3,2)-Mt*Cp*Te8i*Pld(L)/L/(L-l) 

1050  Et(3,l)=Et(3,l)-Mr*Te7r*St*Pld(L)/R0 
1060  Et(3 ,2)=Et(3 ,2)-Mr*T e7i*  S  t*Pld(L)/R0 
1 070  Bt(3, 1  )=Bt(3, 1  )+Mt*Sp*Tb5r*(Pl(L)-Ct*Pld(L)/L/(L- 1 )) 
1080  Bt(3,2)=Bt(3,2)+Mt*Sp*Tb5i*(Pl(L)-Ct*Pld(L)/L/(L-l)) 
1 090  Bt(3 ,1)-Bt(3,l )+Mt*Sp*Tb4r*Pld(L)/L/(L- 1 ) 

1 100  Bt(3,2)=Bt(3,2)+Mt*Sp*Tb4i*Pld(L)/17(L-l) 

1110  Ep(  1 , 1  )=Ep(  1 , 1  )+(2*L- 1 )  *Et(  1,1) 

1 120  Ep(l,2)=Ep(U)+(2*L-l)*Et(l,2) 

1130  Ep(2, 1  )=Ep(2, 1  )+(2*L- 1  )*Et(2, 1 ) 

1 140  Ep(2,2)=Ep(2,2)+(2*L-l)*Et(2,2) 

1 150  Ep(3, 1  )=Ep(3, 1  )+(2*L- 1  )*Et(3, 1 ) 

1 160  Ep(3,2)=Ep(3,2)+(2*L-l)*Et(3^) 

1 170  Bp(l,l)=Bp(l,l)+(2*L-l)*Bt(l,l) 

1180  Bp(  1 ,2)=Bp(  1 ,2)+(2*L- 1  )*Bt(  1 ,2) 

1 190  Bp(2, 1  )=Bp(2, 1  )+(2*L- 1  )*Bt(2, 1 ) 

1200  Bp(2,2)=Bp(2,2)+(2*L- 1  )*B t(2,2) 

1210  Bp(3,l )=Bp(3, 1  )+(2*L- 1  )*Bt(3, 1 ) 

1220  Bp(3,2)=Bp(3,2)+(2*L-l)*Bt(3,2) 

1230  Ner  1  =Et(  1 , 1  )*Et(  1 , 1  )+Et(  1 ,2)*Et(  1 ,2) 

1 240  Ner2=Et(2, 1  )*Et(2, 1  )+Et(2,2)*Et(2,2) 

1 250  Ner3=Et(3, 1  )*Et(3, 1  )+Et(3,2)*Et(3,2) 

1 260  Dne  1  =Ep(  1,1)  *Ep(  1 , 1  )+Ep(  1  ^)*Ep(  1 2) 

1270  Dne2=Ep(2,l)*Ep(2,l)+Ep(2,2)*Ep(2,2) 

1 280  Dne3=Ep(3, 1  )*Ep(3, 1  )+Ep(3,2)*Ep(3,2) 

1290  Nbr  1  =Bt(  1 , 1 )  *  Bt(  1 , 1  )+Bt(  1 ,2)*Bt(  1 ,2) 

1 300  Nbr2=Bt(2, 1  )*Bt(2, 1  )+B t(2,2)*B t(2,2) 

1310  Nbr3=Bt(3,l)*Bt(3,l)+Bt(3,2)*Bt(3,2) 

1 320  Dnb  1  =Bp(  1 , 1  )*  Bp(  1 , 1  )+Bp(  1 ,2)*Bp(  1 ,2) 

1330  Dnb2=Bp(2, 1  )*Bp(2, 1  )+Bp(2,2)*Bp(2,2) 

1 340  Dnb3=Bp(3, 1  )*Bp(3, 1  )+Bp(3,2)*Bp(3,2) 

1350  IF  L=Lmax  THEN  1570 
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1360  IF  Dnbl=0  THEN  1380 

1370  IF  Nbrl/Dnbl<Err*Err  THEN 

1380  IF  Dnb2=0  THEN  1400 

1 390  IF  Nbr2/Dnb2<Err*Err  THEN 

1400  IF  Dnb3=0  THEN  1420 

1410  IF  Nbr3/Dr_b?<Err*Err  THEN 

1420  IF  Dnel=0  THEN  1440 

1430  IFNerl/Dnel<Err*Err  THEN 

1440  IF  Dne2=0  THEN  1460 

1450  IF  Ner2/Dne2<Err*Err  THEN 

1460  IF  Dne3=0  THEN  1480 

1470  IF  Ner3/Dne3<Err*Err  THEN 

1480  GOTO  1570 

1490  END  IF 

1500  END  IF 

1510  END  IF 

1520  END  IF 

1530  END  IF 

1540  END  IF 

1550  L=L+1 

1560  GOTO  550 

1570  Tmpr=Mefacr*Ep(  1 , 1  )-Mefaci*Ep(  1 ,2) 

1 5 80  T mpi=Mefacr*Ep(  1 ,2)+Mefaci*Ep(l ,  1 ) 

1590  Ep(l,l)=-Tmpr 
1600  Ep(l,2)=-Tmpi 

1610  Tmpr=Mmfacr*Bp(l,l)-Mmfaci*Bp(l,2) 

1620  Tmpi=Mmfacr*Bp(l,2)+Mnifaci*Bp(l,l) 

1630  Bp(l,l)=Tmpr 
1640  Bp(U)=TmDi 

1650  Tmpr=Mefacr*Ep(2, 1  )-Mefaci*Ep(2,2) 

1 660  T mpi=Mefacr*Ep(2,2)+Mefaci*Ep(2, 1 ) 

1670  Ep(2,l)=-Tmpr 
1680  Ep(2,2)=-Tmpi 

1 690  Tmpr==Mmfacr*Bp(2, 1  )-Mmfaci*Bp(2,2) 

1700  Tmpi=Mmfacr*Bp(2,2)+Mmfaci*Bp(2,l) 

1710  Bp(2,l)=Tmpr 
1720  Bp(2,2)=Tmpi 

1730  Tmpr=Mefacr*Ep(3,l)-Mefaci*Ep(3,2) 

1740  Tmpi=Mefacr*Ep(3,2)+Mefaci*Ep(3,l) 

1750  Ep(3,l)=Tmpr 
1760  Ep(3,2)=Tmpi 

1770  Tmpr=Mmfacr*Bp(3,l)-Mmfaci*Bp(3,2) 

1780  Tmpi=Mmfacr*Bp(3,2)+Mmfaci*Bp(3,l) 

1790  Bp(3,l)=-Tmpr 
1800  Bp(3^)=-Tmpi 

1810  CALL  Geomfldb_e(Q,St,Cp,Sp,Bp(*)3p(*)^ds(*),B(*),E(*)) 
1820  PRINT  E(*)  lElectric  field  in  original  frame 
1830  PRINT 

1840  PRINT  B(*)  ’.Magnetic  field  in  original  frame 
1850  END 
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PROGRAM  CURRDIPSPH 
10  OPTION  BASE  1 

20  DIM  Xf(3),Xd(3),P(3),Rds(3,3),Bp(3,2),Ep(3,2))B(3,2)>E(3,2) 

30  DIM  Pl(100)d5ld(100),Ndma(100,2),Ndea(100,2)3t(3,2),Bt(3,2) 

40  INPUT  "FREQUENCY  AND  MEDIUM  AND  SPHERE(Ss<0  =>  INF) 

CONDUCnvmES(MHO/M)?",F,Sm,Ss 

50  INPUT  "RELATIVE  PERMEABILITY  OF  MEDIUM  AND  SPHERE?" ,Mm,Ms 
60  INPUT  "RELATIVE  PERMITTIVITY  OF  MEDIUM  AND  SPHERE?"  JEm,Es 
70  INPUT  "POSITION  OF  FIELD  POINT? (M)",Xf(*) 

80  INPUT  "POSITION  OF  SOURCE  POINT?(M)",Xd(*) 

90  INPUT  "CURRENT  DIPOLE  MOMENT  VECTOR?(AMP-M)",P(*) 

100  INPUT "  SPHERE  RADIUS(M)  AND  MAXIMUM  POLAR  ANGLE  INDEX?",AJE11 

110  Err=l.E-6 

120  Lmax=Ell+l 

130  W=2*PI*F 

140  M0=4*PI*l.E-7 

150  E0=8.85415E-12 

160  Utl=W*M0 

170  Ut2=W*E0*Utl 

180  K12r=Ut2*Es*Ms 

190  K12i=Utl*Ms*Ss 

200  K22r=Ut2*Em*Mm 

210  K22i=Utl*Mm*Sm 

220  Md=SQR(K12r*K12r+K12i*K12i) 

230  K 1  r=SQR((Md+K  1 2r)/2) 

240  Kli=SQR(ABS(Md-K12r)/2) 

250  Md=SQR(K22r*K22r+K22i*K22i) 

260  K2r=SQR((Md+K22r)/2) 

270  K2i=SQR( AB  S  (Md-K22r)/2) 

280  Mk2=K2r*K2r+K2i*K2i 
290  Ik2r=K2r/Mk2 
300  Ik2i=-K2i/Mk2 

310  Eefac=M0*Mm*W/4/PI  !  Electric  fields  in  volt/meter 

320  Eefacr=Eefac*Ik2r 
330  Eefaci=Eefac*Ik2i 

340  Emfac=M0*Min/4/PI  !  Magnetic  fields  in  tesla 

350  Emfacr=-Emfac*K2i 
360  Emfaci=Emfac*K2r 
370  Mur=Ms/Mm 

380  .'Only  frequency  and  medium  constants  to  here 
390  REDIM  Ndma(Lmax,2),Ndea(Lmax,2) 

400  TF  <sc>0  THFN  440 

410  CALL  Ndmicon(Lmax,K2r,K2i,A,Ndma(*)) 

420  CALL  Ndeicon(Lmax,K2r,K2i,A,Ndea(*)) 

430  GOTO  460 

440  CALL  Ndm(Lmax,K  1  r,K  1  i,K2r,K2i,Mur,A,Ndma(*)) 

450  CALL  Nde(Lmax,Klr,Kli,K2r,K2i,Mur,A,Ndea(*)) 

460  !  Added  only  sphere  radius  to  here 

470  CALL  Geomdipsph(Xd(*)T>(*),Rds(*),R0d>r3) 

480  CALL  Geomfldpos(Xf(*),Rds(*),R,Ct,St,Cp,Sp) 

490  Zr=K2r*R 
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500  Zi=K2i*R 

510  Z0r^K2r*R0 

520  Z0i=K2i*R0 

530  REDIM  Pl(Lmax)JPld(Lmax) 

540  CALL  Pl_pldot(Lmax- 1  ,Q,Pl(*)rPld(*)) 

550  MAT  Ep=  (0) 

560  MAT  Bp=  (0) 

570  L=2 

580  MAT  Et=  (0) 

590  MAT  Bt=  (0) 

600  CALL  Hcomb(L-l,Zr,Zi,RrHzr,Hzi,Chzr,Chzi) 

610  CALL  Hcomb(L- 1  ,ZOr,ZOi,RO,HzOr>HzOi>ChzOr,ChzOi) 
620  Tblr=Hzr*HzOr-Hzi*HzOi 
630  Tbli=Hzr*HzOi+Hzi*HzOr 
640  Tb2r=Tb  lr*Ndea(L,  1  )-Tb  1  i*Ndea(L,2) 

650  Tb2i=Tb  lr*Ndea(L,2)+Tb  li*Ndea(L,  1 ) 

660  Te2r=Tb  lr*Ndma(L,  1  )-Tb  1  i*Ndma(L,2) 

670  Te2i=Tblr*Ndma(L,2)+Tbli*Ndma(L,l) 

680  Bt(l,l)=Pt*St*Sp*Tb2r*Pld(L)/R 
690  Bt(U)=Pt*St!Sp*Tb2i*Pld(L)/R 
700  Tb3r=Hz0r*Chzr-Hz0i*Chzi 
710  Tb3i=Hz0r*Chzi+Hz0i*Chzr 
720  Tb4r=Tb3r*Ndea(L,  1  )-Tb3i*Ndea(L,2) 

730  Tb4i=Tb3r*Ndea(L,2)+Tb3i*Ndea(L,l) 

740  Tb5r=Hzr*ChzOr-Hzi*ChzOi 
750  Tb5i=Hzr*Chz0i+Hzi*Chz0r 
760  Tb6r=Tb5r*Ndma(L,l)-Tb5i*Ndma(L^2) 

770  Tb6i=Tb5r*Ndma(L,2)+Tb5i*Ndma(L,  1) 

780  Et(  1 , 1  )=Pr*Te2r*L*(L-  1)*P1(L)/R/R0 
790  Et(  1 ,2)=Pr*Te2i*L*(L-  1)*P1(L)/R/R0 
800  Et(  1,1  )=Et(  1 , 1  )+Pt*Cp*St*Tb6r*Pld(L)/R 
8 10  Et(  1 ,2)=Et(  1 ,2)+Pt*Cp*St*Tboi*Pld(L)/R 
820  Te3r=Chzr*Chz0r-Chzi*Chz0i 
830  Te3i=Chzr*Chz0i+Chzi*ChzQr 
840  Te4r=Te3r*Ndma(L,  l)-Te3i*Ndma(L,2) 

850  Te4i=Te3r*Ndma(L,2)+Te3i*Ndma(L,  1) 

860  Te5r==-K22r*Tb2r+K22i*Tb2i 

870  Te5i=-K22r*Tb2i-K22i*Tb2r 

880  Te6r=Tb3r*Ndma(L,  l)-Tb3i*Ndma(L^) 

890  Te6i=Tb3r*Ndma(L,2)+Tb3i*Ndma(L,l) 

900  Bt(2, 1  )=Bt(2, 1  )+Pt*Sp*Tb4r*(Pl(L)-Ct*Pld(L)/IV(L- 1 )) 
910  Bt(2,2)=Bt(2,2)+Pt*Sp*Tb4i*(Pl(L)-Ct*Pld(L)/L/(L-l)) 
920  Bt(2,l)=Bt(2, 1  )-Pt*Sp*Tb6r*Pld(L)/L/(L- 1 ) 

930  Bt(2,2)=Bt(2,2)-Pt*Sp*Tb6i*Pld(L)/L/(L- 1 ) 

940  Et(2,l  )=Et(2, 1  )+Pr*S  t*Te6r*Pld(L)/R0 
950  Et(2,2)=Et(2,2)+Pr*St*Te6i*Pld(L)/R0 
960  Et(2,l)=Et(2,l)+Pt*Cp*Te5r*Pld(L)/L/(L-l) 

970  Et(2,2)=Et(2,2)+Pt*Cp*Te5i*Pld(L)/L/(L-l) 

980  Et(2, 1  )=Et(2, 1  )-Pt*Cp*Te4r*(Pl(L)-Q*Pld(L)/L/(L-l» 
990  Et(2,2)=Et(2,2)-Pt*Cp*Te4i*(Pl(L)-Ct*Pld(L)/L/(L-l)) 
1000  Tb7r=Tblr*Ndma(L,l)-Tbli*Ndma(L^) 

1010  Tb7i=Tblr*Ndma(L,2)+Tbli*Ndma(L,l) 
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1020  Tb8r=Tb3r*Ndea(L,  1  )-Tb3i*Ndea(L,2) 

1030  Tb8i=Tb3r*Ndea(L,2)+Tb3i*Ndea(L,l) 

1040  Bt(3,l)=Bt(3,l)-Pt*Cp*Tb6r*(Pl(L)-Ct*Pld(L)/L/(L-l)) 
1050  Bt(3,2)=Bt(3,2)-Pt*Cp*Tb6i*(Pl(L)-Ct*Pld(L)/L/(L-l)) 
1060  Bt(3,l)=Bt(3,l)+Pt*Cp*Tb8r*Pld(L)/L/(L-l) 

1070  Bt(3 ,2)=B t(3,2)+Pt*Cp*Tb8i*Pld(L)/L/(L- 1 ) 

1080  Bt(3, 1 ) =B t(3 , 1  )+Pr*Tb7 r * S t* Pld(L)/R0 
1090  Bt(3,2)=Bt(3,2)+Pr*Tb7i*St*Pld(L)/R0 
1 100  Et(3 , 1  )=Et(3, 1  )-Pt*Sp*Te5r* (Pl(L)-Ct*Pld(L)/LV (L- 1 )) 
1 1 10  Et(3,2)=Et(3,2)-Pt*Sp*Te5i*(Pl(L)-Ct*Pld(L)/U(L-l)) 
1 120  Et(3,l)=Et(3,l)+Pt*Sp*Te4r*Pld(L)/L/(L-l) 

1 130  Et(3,2)=Et(3,2)+Pt*Sp*Te4i*Pld(L)AV(L-l) 

1 1 40  Bp(  1 , 1  )=Bp(  1 , 1  )+(2*L- 1  )*Bt(  1 , 1 ) 

1 150  Bp(  1 ,2)=Bp(  1 ,2)+(2*L- 1  )*Bt(  1 ,2) 

1 1 60  Bp(2, 1  )=Bp(2, 1  )+(2*L- 1  )*B t(2, 1 ) 

1 170  Bp(2,2)=Bp(2,2)+(2*L-l)*Bt(2,2) 

1 1 80  Bp(3 , 1  )=Bp(3, 1  )+(2*L- 1  )*Bt(3, 1 ) 

1 190  Bp(3,2)=Bp(3,2)+(2*L-l)*Bt(3,2) 

1200  Ep(l,l)=Ep(l,l)-(2*L-l)*Et(l,l) 

1 2 1 0  Ep(  1 ,2)=Ep(  1 ,2)-(2*L- 1  )*Et(  1 ,2) 

1 220  Ep(2, 1  )=Ep(2, 1  )+(2*L- 1 )  *Et(2, 1 ) 

1230  Ep(2,2)=Ep(2,2)+(2*L-l)*Et(2,2) 

1 240  Ep(3, 1  )=Ep(3, 1  )+(2*L- 1  )*Et(3, 1 ) 

1250  Ep(3,2)=Ep(3,2)+(2*L-  l)*Et(3,2) 

1260  Ner  1  =Et(  1 , 1  )*Et(  1 , 1  )+Et(  1 ,2)*Et(  1 ,2) 

1 270  Ner2-Et(2, 1  )*Et(2, 1  )+Et(2,2)*Et(2,2) 

1280  Ner3=Et(3 , 1  )*Et(3 , 1  )+Et(3 ,2)  *Et(3,2) 

1290  Dnel=Ep(l,l)*Ep(l,l)+Ep(l,2)*Ep(U) 

1300  Dne2=Ep(2,l)*Ep(2,l)+Ep(2^)*Ep(2^) 

1310  Dne3=Ep(3,l)*Ep(3,l)+Ep(3,2)*Ep(3,2) 

1 320  Nbr  1  =Bt(  1 , 1  )*Bt(  1 , 1  )+B t(  1 ,2)*Bt(  1 ,2) 

1330  Nbr2=Bt(2,l)*Bt(2,l)+Bt(2,2)*Bt(2,2) 

1340  Nbr3=Bt(3,l)*Bt(3,l)+Bt(3,2)*Bt(3,2) 

1350  Dnb  1  =Bp(  1 , 1  )*Bp(  1 , 1  )+Bp(  1  ^)*Bp(  1 ,2) 

1360  Dnb2=Bp(2,l)*Bp(2,l)+Bp(2,2)*Bp(2^2) 

1370  Dnb3=Bp(3,l)*Bp(3,l)+Bp(3,2)*Bp(3,2) 

1380  IF  L=Lmax  THEN  1600 
1390  IFDnbl=0THEN  1410 
1400  IF  Nbrl/Dnbl<Err*Err  THEN 
1410  IF  Dnb2=0  THEN  1430 
1420  IF  Nbr2/Dnb2<Err*Err  THEN 
1430  IF  Dnb3=0  THEN  1450 
1440  EF  Nbr3/Dnb3<Err*Err  THEN 
1450  IF  Dnel=0  THEN  1470 

1460  IF  Nerl/Dnel<Err*Err  THEN 

1470  IF  Dne2=0  THEN  1490 

1480  EF  Ner2/Dne2<Err*Err  THEN 

1490  IF  Dne3=0  THEN  1510 

1500  IF  Ner3/Dne3<Err*Err  THEN 

1510  GOTO  1600 
1520  END  IF 
1530  END  IF 
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1540  END  IF 
1550  END  IF 
1560  END  IF 
1570  END  IF 
1580  L=L+1 
1590  GOTO  580 

1600  Tmpr=Emfacr*Bp(  1 , 1  )-Emfaci*Bp(  1 ,2) 

1610  Tmpi=Emfacr*Bp(  1 ,2)+Emfaci*Bp(  1,1) 

1620  Bp(l,l)=Tmpr 
1630  Bp(l,2)=Tmpi 

1 640  Tmpr=Eefacr*Ep(  1 , 1  )-Eefaci*Ep(  1 ,2) 

1 650  Tmpi=Eefacr*Ep(l  ,2)+Eefaci*Ep(  1 , 1 ) 

1660  Ep(l,l)=Tmpr 
1670  Ep(l,2)=Tmpi 

1680  Tmpr=Emfacr*Bp(2, 1  )-Emfaci*Bp(2,2) 

1690  Tmpi=Emfacr*Bp(2,2)+Emfaci*Bp(2,l) 

1700  Bp(2,l)=Tmpr 
1710  Bp(2,2)=Tmpi 

1720  Tmpr=Eefacr*Ep(2,l)-Eefaci*Ep(2,2) 

1730  Tmpi=Eefacr*Ep(2,2)+Eefaci*Ep(2,l) 

1740  Ep(2,l)=Tmpr 
1750  Ep(2,2)=Tmpi 

1760  Tmpr=Emfacr*Bp(3,l)-Emfaci*Bp(3,2) 

1770  Tmpi=Emfacr*Bp(3,2)+Emfaci*Bp(3, 1) 

1780  Bp(3,l)=Tmpr 
1790  Bp(3,2)=Tmpi 

1 800  Tmpr=Eefacr*Ep(3,l)-Eefaci*Ep(3,2) 

1810  Tmpi=Eefacr*Ep(3,2)+Eefaci*Ep(3, 1 ) 

1820  Ep(3,l)=Tmpr 
1830  Ep(3,2)=Tmpi 

1840  CALL  Geomfldb_e(Ct,St,Cp,Sp,Bp(*),Ep(*),Rds(*),B(*),E(*)) 
1850  PRINT  E(*)  ! Electric  field  in  original  frame 
1860  PRINT 

1870  PRINT  B(*)  ! Magnetic  field  in  original  frame 
1880  END 


PROGRAM  MGDIPSPHDC 
10  OPTION  BASE  1 

20  DIM  Xf(3),Xd(3),M(3),Rds(3,3),Bp(3),Ep(3),B(3),E(3) 

30  DEMPl(100)d>ld(100)3t(3),Bt(3) 

40  INPUT  "MEDIUM  AND  SPHERE  CONDUdTVmES(MHO/M)?",Sm,Ss 
50  INPUT  "RELATIVE  PERMEABILITY  OF  MEDIUM  AND  SPHERE?\Mm,Ms 
60  INPUT  "POSITION  OF  FIELD  POINT?(M)",Xf(*) 

70  INPUT  "POSITION  OF  SOURCE  POINT?(M)"^d(*) 

80  INPUT  "MAGNETIC  DIPOLE  MOMENT  VECTOR?(AMP-MA2)",M(*) 

90  INPUT "  SPHERE  RADIUS(M)  AND  MAXIMUM  POLAR  ANGLE  INDEX?", AJEU 

100  Err=l.E-6 

110  Lmax=Ell+l 

120  T=Ms/Mm 

130  M0=4*PI*l.E-7 
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140  !Only  medium  constants  to  here 

150  CALL  Geomdipsph(Xd(*),M(*),Rds(*),RO,Mr,Mt) 

160  CALL  Geomfldpos(Xf(*),Rds(*),R,Ct,St,Cp,Sp) 

170  Facm=M0*Mm*(T-l)/4/PI/A/R/R0 
180  REDIM  Pl(Lmax),Pld(Lmax) 

190  CALL  Pl_pldot(Lmax- 1  ,Ct,Pl(*),Pld(*)) 

200  MAT  Ep=  (0) 

210  MAT  Bp=  (0) 

220  Rat=A*A/R/R0 
230  Ratl=Rat 
240  L=2 
250  MAT  Bt=  (0) 

260  Ratl=Ratl*Rat 
270  Dt=l/((T+1)*(L-1)+1) 

280  Bt(l)=-Mt*St*Cp*L*(L-l)*Ratl*Pld(L)*Dt 

290  Bt(l)=Bt(l)+Mr*L*L*(L-l)*Ratl*Pl(L)*Dt 

300  Bt(2)=Mt*Cp*(L-l)*Ratl*(L*(L-l)*Pl(L)-Ct*Pld(L))*Dt 

310  Bt(2)=Bt(2)+Mr*St*L*(L-l)*Pld(L)*Ratl*Dt 

320  Bt(3)=-Mt*Sp*(L-l)*Rad*Pld(L)*Dt 

330  MAT  Bp=  Bp+Bt 

340  Nbrl=Bt(l)*Bt(l) 

350  Nbr2=Bt(2)*Bt(2) 

360  Nbr3=Bt(3)*Bt(3) 

370  Dnb  1  =Bp(l)*Bp(  1) 

380  Dnb2=Bp(2)*Bp(2) 

390  Dnb3=Bp(3)*Bp(3) 

400  IF  L=Lmax  THEN  530 

410  IF  Dnb  1=0  THEN  430 

420  IF  Nbrl/Dnbl<Err*Err  THEN 

430  IF  Dnb2=0  THEN  450 

440  IF  Nbr2/Dnb2<Err*Err  THEN 

450  IF  Dnb3=0  THEN  470 

460  IF  Nbr3/Dnb3<Err*Err  THEN 

470  GOTO  530 

480  END  IF 

490  END  IF 

500  END  IF 

510  L=L+1 

520  GOTO  250 

530  Bp(  1  )=Facm*Bp(  1 ) 

540  Bp(2)=Facm*Bp(2) 

550  Bp(3)=Facm*Bp(3) 

560  CALL  Geomfldb_e(a,St,Cp,Sp,Bp(*)3p(*)^ds(*),B(*)3(*)) 
570  PRINT  E(*)  !DC  electric  field  in  original  frame(volt/meter) 
580  PRINT 

590  PRINT  B(*)  !DC  magnetic  field  in  original  frame(tesla) 

600  END 
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PROGRAM  CRDIPSPHDC 
10  OPTION  BASE  1 

20  DIM  Xf(3),Xd(3),P(3),Rds(3,3),Bp(3),Ep(3),B(3),E(3) 

30  DIM  Pl(100),Pld(100),Et(3),Bt(3) 

40  INPUT  "MEDIUM  AND  SPHERE  CONDUCnvnTES(MHO/M)?'\Sm,Ss 
50  INPUT  "RELATIVE  PERMEABILITY  OF  MEDIUM  AND  SPHERE?", Mm,Ms 
60  INPUT  "POSITION  OF  FIELD  POINT?(M)",Xf(*) 

70  INPUT  "POSITION  OF  SOURCE  POINT?(M)"yXd(*) 

80  INPUT  "CURRENT  DIPOLE  MOMENT  VECTOR?(AMP-M)",P(*) 

90  INPUT "  SPHERE  RADIUS(M)  AND  MAXIMUM  POLAR  ANGLE  INDEX?", AJE11 

100  Err=l.E-6 

110  Lmax=Ell+l 

120  M0=4*PI*l.E-7 

130  Facm=M0*Mm/4/PI/A 

140  Face=l/4/PI/A/Sm 

141  T= Ms/Mm 
150  Ee=Ss/Sm 

170  !Only  medium  constants  to  here 

180  CALL  Geomdipsph(Xd(*)J>(*),Rds(*),ROJPrJ>t) 

190  CALL  Geomfldpos(Xf(*),Rds(*),R,Ct,St,Cp,Sp) 

200  REDIM  Pl(Lmax)JPld(Lmax) 

210  CALL  Pl_pldot(Lmax-l,Ct,Pl(*)JPld(*)) 

220  MAT  Ep=  (0) 

230  MAT  Bp=  (0) 

240  Rat=A*A/R/R0 
250  Ratl=Rat 
260  L=2 
270  MAT  Et=  (0) 

280  MAT  Bt=  (0) 

290  Ratl=Ratl*Rat 
300  Dt=l/((T+1)*(L-1)+1) 

310  De=l/((Ee+l)*(H)+l) 

320  Bt(  1  )=(T- 1  )*Pt*St*Sp*L*Ratl*Pld(L)*Dt 

330  Et( l)=(Ee-  l)*Pr*L*L*(L- l)*Pl(L)*Ratl*De 

340  Et(l)=Et(l)-(Ee-l)*Pt*Cp*St*L*(L-l)*Ratl*Pld(L)*De 

350  Bt(2)=Pt*Sp*(Ee-l)*Ratl*Pld(L)/R0*De 

360  Bt(2)=Bt(2)-Pt*Sp*(T-l)*(L*(L-l)*Pl(L)-Ct*Pld(L))*Ratl/R*Dt 

370  Et(2)=(Ee- 1  )*Pr*S  t*L*(L- 1  )*Ratl  *Pld(L)*De 

380  Et(2)=Et(2)+(Ee- 1 )  *Pt*Cp*(L- 1 )  *Ratl*(L*  (L- 1  )*Pl(L)-Ct*Pld(L))*De 

390  Bt(3)=(Ee-l)*Pr*St*L*Ratl*Pld(L)*De/R0 

400  Bt(3)=Bt(3)+(Ee- l)*Pt*Cp*(L*(L-l)*Pl(L)-Q*Pld(L))*Ratl/RO*De 

410  Bt(3)=Bt(3)-(T-l)*Pt*Cp*Ratl*Pld(L)/R*Dt 

420  Et(3)=-Pt*Sp*(Ee-  1)*(L- 1  )*Ratl*Pld(L)*De 

430  MAT  Bp=  Bp+Bt 

440  MAT  Ep=  Ep+Et 

450  Nerl=Et(l)*Et(l) 

460  Ner2=Et(2)*Et(2) 

470  Ner3=Et(3)*Et(3) 

480  Dnel=Ep(l)*Ep(l) 

490  Dne2=Ep(2)*Ep(2) 

500  Dne3=Ep(3  )*Ep(3) 
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510  Nbrl=Bt(l)*Bt(l) 

520  Nbr2=Bt(2)*Bt(2) 

530  Nbr3=Bt(3)*Bt(3) 

540  Dnb  1  =Bp(  1  )*Bp(  1 ) 

550  Dnb2=Bp(2)*Bp(2) 

560  Dnb3=Bp(3)*Bp(3) 

570  IF  L=Lmax  THEN  790 

580  IF  Dnb  1=0  THEN  600 

590  IF  Nbrl/Dnbl<Err*Err  THEN 

600  IF  Dnb2=0  THEN  620 

610  IF  Nbr2/Dnb2<Err*Err  THEN 

620  IF  Dnb3=0  THEN  640 

630  IF  Nbr3/Dnb3<Err*Err  THEN 

640  IF  Dnel=0  THEN  660 

650  IF  Nerl/Dnel<Err*Eir  THEN 

660  IF  Dne2=0  THEN  680 

670  IF  Ner2/Dne2<Err*Eir  THEN 

680  IF  Dne3=0  THEN  700 

690  IF  Ner3/Dne3<Err*Err  THEN 

700  GOTO  790 

710  END  IF 

720  END  IF 

730  END  IF 

740  END  IF 

750  END  IF 

760  END  IF 

770  L=L+1 

780  GOTO  270 

790  Bp(  1  )=Facm*Bp(  1)/R 

800  Ep(  1  )=Face*Ep(  1  )/R/R0 

810  Bp(2)=Facm*Bp(2) 

820  Ep(2)=Face*Ep(2)/R/R0 
830  Bp(3)=Facm*Bp(3) 

840  Ep(3)=Face*Ep(3)/R0/R 

850  CALL  Geomfldb_e(Ct,St,C,  ,Sp3?(*),Ep(*),Rds(*),B(*)JE(*)) 
860  PRINT  E(*)  !DC  electric  field  in  original  frame  (volt/meter) 
870  PRINT 

880  PRINT  B(*)  !DC  magnetic  field  in  original  firame(tesla) 

890  END 


10  SUB  Geomdipsph(D(*),V(*)3(*)»Din,Vr,Vt) 

20  ! 

30  ! Produces  transformation  R(*)  from  general  frame  with  source  com- 

40  Iponents  D(*)  and  moment  components  V(*)  to  the  standard  frame 
50  !with  source  on  the  z-axis,  radial  component  Vr,  and  tangential 
60  Icomponent  Vt  oriented  along  the  positive  x-axis. 

70  ! 

80  OPTION  BASE  1 

90  DIM  Rl(3,3),R2(3,3),R3(3,3),Rt(3,3),Vtm(3) 

100  MAT  Rl=  (0) 

110  MAT  R2=  (0) 
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120  MAT  R3=  (0) 

130  Dh=SQR(D(  1  )*D(  1)+D(2)*D(2)) 

140  Dm=SQR(Dh*Dh+D(3)*D(3)) 

150  Std=Dh/Dm 
160  Ctd=D(3)/Dm 
170  IF  Std>l.E-12  THEN 
180  Cpd=D(l)/Dh 
190  Spd=D(2)/Dh 
200  Rl(l,l)=Cpd 
210  Rl(U)=Spd 
220  Rl(2,l)=-Spd 
230  Rl(2,2)=Cpd 

240  Rl(3,3)=l 

250  R2(l,l)=Ctd 
260  R2(l,3)=-Std 

270  R2(2,2)=l 

280  R2(3,l)=Std 

290  R2(3,3)=Ctd 

300  MAT  Rt=  R2*R1 
310  ELSE 
320  MAT  Rt=  IDN 
330  END  IF 
340  MAT  Vtm=  Rt*  V 

350  Vh=SQR(Vtm(l)*Vtm(l)+Vtm(2)*Vtm(2)) 

360  Vm=SQR(Vh*Vh+Vtm(3)*Vtm(3)) 

370  IF  Vh/Vm>l.E-12  THEN 

380  Cpv=Vtm(l)/Vh 

390  Spv=Vtm(2)/Vh 

400  R3(l,l)=Cpv 

410  R3(l,2)=Spv 

420  R3(3,l)=0 

430  R3(2,l)=-Spv 

440  R3(2^2)=Cpv 

450  R3(3,3)=l 

460  ELSE 

470  MAT  R3=  IDN 

480  END  IF 

490  MAT  R=  R3*Rt 

500  Vt=Vh 

510  Vr=Vtm(3) 

520  SUBEND 

10  SUB  Geomfldpos(Xf(*),Rds(*),R,Ct,St,Cp,Sp) 

20  ! 

30  ’.Transforms  field  point  in  general  frame  to  special  dipole- 
40  !  Sphere  frame  and  converts  it  to  polar  coordinates  in  that 
50  .’Frame. 

60  ! 

70  OPTION  BASE  1 
80  DIM  X(3) 

90  MAT  X=  Rds*Xf 

100  Rh=SQR(X(  1  )*X(  1  )+X(2)*X(2)) 
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110  R=SQR(Rh*Rh+X(3)*X(3)) 

120  St=Rh/R 

130  Q=X(3)/R 

140  IF  St>l-E-12  THEN 

150  Cp=X(l)/Rh 

160  Sp=X(2)/Rh 

170  ELSE 

180  Cp=l 

190  Sp=0 

200  END  IF 

210  SUBEND 

10  SUB  Geomfldb_e(Ct,St,Cp^p3p(*)»Ep(*)>Rds(*)3(*)»E(*)) 
20  ! 

30  {Transforms  complex  polar  coordinate  fields  in  the  special 
40  IDipole-sphere  frame  to  cartesian  coordinates  in  that  frame 
50  !  And  then  rotates  the  fields  back  to  the  original  general 

60  {Frame. 

70  ! 

80  OPTION  BASE  1 

90  DIM  Irds(3,3)TPc(3,3),Et(3,2)3t(3t2) 

100  MAT  Irds=  INV(Rds) 

110  Pc(l.D=St*Cp 
120  Pc(U)=Q*Cp 
130  Pc(l,3)=-Sp 

140  Pc(2,l)=St*Sp 
150  Pc(2,2)=Ct*Sp 

160  Pc(2,3)=Cp 
170  Pc(3,l)=Ct 
180  Pc(3,2)=-St 

190  Pc(3,3)=0 

200  MAT  Bt=  Pc*Bp 
210  MAT  Et=  Pc*Ep 
220  MAT  B=  Irds*Bt 
230  MAT  E=  Irds*Et 
240  SUBEND 

10  SUB  Jcomb(L,Ur,UiJUr JiJcrJd) 

20  CALL  Spherejnz(L,Ur,UiJrJi) 

30  CALL  Spherejnz(L+ 1  ,Ur,Ui,J  lrj  1  i) 

40  Jcr=((L+l)*  Jr-(Ur*  J  lr-Ui*Jli))/R 
50  Jci=((L+ 1  )*  Ji-(Ur*  J 1  i+Ui*  J  lr))/R 
60  SUBEND 

10  SUB  Hcomb(L,Ur,Ui3,Hr4B4Icr,Hci) 

20  CALL  Spherehnz(L,Ur,Ui3r>Hi) 

30  CALL  Spherehnz(L+l,Ur,Ui,Hlr>Hli) 

40  Hcr=((L+l)*Hr-(Ur*Hlr-Ui*Hli))/R 
50  Hci=((L+l)*Hi-(Ur*Hli+Ui*Hlr))/R 
60  SUBEND 
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10  SUB  Spherejnz(N,Zr,Zi,Sjr,Sji) 

20  CALL  Jnuevrywhr(N+.5,Zr,Zi,JrJi) 

30  Dm=Zr*Zr+Zi*Zi 
40  Ur=Zr/Dm 
50  Ui=-Zi/Dm 
60  Nm=SQR(Ur*Ur+Ui*Ui) 

70  Sur=SQR((Nm+Ur)/2) 

80  Sui=SGN(Ui)*SQR(ABS(Nm-Ur)/2) 

90  Sjr=SQR(PI/2)*(Sur*Jr-Sui*Ji) 

100  Sji=SQR(PI/2)*(Sur*Ji+Sui*Jr) 

110  SUBEND 

10  SUB  Spherehnz(N,Zr,Zi,Shr,Shi) 

20  CALL  Hnuevry  whr(N+.5  ,Zr,Zi  Jir  Jli) 

30  Dm=Zr*Zr+Zi*Zi 
40  Ur=Zr/Dm 
50  Ui=-Zi/Dm 
60  Nm=SQR(Ur*Ur+Ui*Ui) 

70  Sur=SQR((Nm+Ur)/2) 

80  Sui=SGN(Ui)*SQR(ABS(Nm-Ur)/2) 

90  Shr=SQR(PI/2)*(Sur*Hr-Sui*Hi) 

100  Shi=SQR(PI/'2)*(Sur*Hi+Sui*Hr) 

110  SUBEND 

10  SUB  Jnuevrywhr(Nu,Zr,Zi,Jr  Ji) 

20  Zmag=SQR(Zr*Zr+Zi*Zi) 

30  IF  Nu<10  THEN 

40  IF  ZmagclO  THEN 

50  CALLJnu(Nu,Zr^iJrJi,l£-28) 

60  ELSE 

70  IF  Nu=0  THEN 

80  CALL  Jnasy(Nu,Zr,Zi,Jr  Ji,5) 

90  ELSE 

100  CALL  Jfiuniasym(Nu,Zr,Zi,Jr,Ji,5) 

110  END  IF 

120  END  IF 
130  ELSE 

140  CALL  Jfiuniasym(Nu,Zr,Zi  Jr,Ji,5) 

150  END  IF 
160  SUBEND 

10  SUB  Hnuevrywhr(Nu,Zr,Zi,Hr,Hi) 

20  Zmag=SQR(Zr*Zr+Zi*Zi) 

30  IFNu<10THEN 

40  IF  ZmagclO  THEN 

50  CALLHlnu(Nu,Zr,Zi,Hr>Hi,l.E-28) 

60  ELSE 

70  IF  Nu=0  THEN 

80  CALL  Hlnasy(Nu,Zr,Zi,Hr>Hi,5) 

90  ELSE 

100  CALL  Hfkiimasym(Nu,Zr,Zi,Hr,Hi,5) 

110  END  IF 
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120  END  IF 
130  ELSE 

140  CALL  Hfkuniasym(Nu,Zr,Zi,Hr,Hi,5) 

150  END  IF 
160  SUBEND 

10  SUB  PI  _pldot(N,X,Pl(*),Pldot(*)) 

20  P1(1)=T 

30  IF  N=0  THEN  60 

40  P1(2)=X 

50  IF  N=1  THEN  80 
60  Pldot(l)=0 
70  IF  N=0  THEN  170 

80  Pldot(2)=l 
90  IFN=1  THEN  170 
100  FOR  L=3  TO  N+l 
110  C1=1/(L-1) 

120  C2=1-C1 

130  C3=2-C1 

140  P1(L)=C3*X*P1(L- 1  )-C2*Pl(L-2) 

150  Pldot(L)=(C3  *X*Pldot(L- 1  )~Pldot(L-2))/C2 

160  NEXT  L 
170  SUBEND 

10  SUB  Nde(Lmax,Klr,Kli,K2r,K2i,Miir,R,Ndea(*)) 

20  Ulr=Klr*R 

30  Uli=KU*R 

40  U2r=K2r*R 

50  U2i=K2i*R 

60  FOR  L=1  TOLmax 

70  ON  ERROR  GOTO  1 10 

80  CALL  Spherejnz(L- 1  ,U  lr,U  li  J  lr,Jli) 

90  CALL  Spherejnz(L,U  lr.UliJ  1  lr,Jl  li) 

100  GOTO  140 

110  Jlratr=0 

120  Jlrati=l 

130  GOTO  190 

140  Dm=Jlr*Jlr+Jli*Jli 

150  Rpr=Jlr/Dm 

160  Rpi=-Jli/Dm 

170  Jlratr=Rpr*Jl  lr-Rpi*Jl  li 

180  Jlrati=Rpr*Jlli+Rpi*Jllr 

190  OFF  ERROR 

200  CALL  Spherejnz(L- 1  ,U2r,U2iJ2rJ2i) 

210  CALLSpherejnz(L,U2r,U2U21rJ21i) 

220  CALL  Spherehnz(L- 1  ,U2r,U2i,H2r,H2i) 

230  CALL  Spherchnz(L,U2r,U2idI21rJI21i) 

240  Nblr=L*J2r-(U2r*J21r-U2i*J21i) 

250  Nbli=L*J2i-(U2r*J21i+U2i*J21r) 

260  Nb2r=L-(Ulr*Jlratr-Uli*Jlrati) 

270  Nb2i=-(U  li +  J  lrati+U  li*  Jlratr) 

280  Topr=Mur*Nblr-(J2r*Nb2r-J2i*Nb2i) 
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290  Topi=Mur*Nbli-(J2r*Nb2i+J2i*Nb2r) 

300  Db  1  r=L*H2r-(U2r*H2  lr-U2i*H2 1  i) 

310  Db  1  i=L*H2i-(U  2r*H2 1  i+U2i*H2  lr) 

320  Btmr=Mur*Dblr-(H2r*Nb2r-H2i*Nb2i) 

330  Btmi=Mur*Dbli-(H2r*Nb2i+H2i*Nb2r) 

340  Bm=Btmr*Btmr+Btmi*Btmi 

350  Rbmr=Btmr/Bm 

360  Rbmi=-Btmi/Bm 

370  Ndea(L,l)=-(Topr*Rbmr-Topi*Rbmi) 

380  Ndea(L,2)=-(Topr*Rbmi+Topi*Rbmr) 

390  NEXT  L 
400  SUBEND 

10  SUB  Ndm(Lmax,Klr,Kli,K2r,K2i,Mur^,Ndma(*)) 

20  Ulr=Klr*R 

30  Uli=Kli*R 

40  U2r=K2r*R 

50  U2i=K2i*R 

60  Klsqr=Klr*Klr-Kli*Kli 

70  Klsqi=2*Klr*Kli 

80  K2sqr=K2r*K2r-K2i*K2i 

90  K2sqi=2*K2r*K2i 

100  Dm=Klr*Klr+Kli*Kli 

110  Rpr=Klr/Dm 

120  Rpi=-Kli/Dm 

130  Ep2r=Rpr*K2r-Rpi*K2i 

140  Ep2i=Rpr*K2i+Rpi*K2r 

150  FOR  L=1  TO  Lmax 

160  ON  ERROR  GOTO  200 

170  CALL  Spherejnz(L- 1  ,U  lr,U  li  J  lr,Jli) 

180  CALL  Spherejnz(L,U  lr,U  li  J1  lr,Jl  li) 

190  GOTO  230 

200  Jlratr=0 

210  Jlrati=l 

220  GOTO  280 

230  Dm=J  lr*J  lr+Jli*Jli 

240  Rpr=Jlr/Dm 

250  Rpi=-Jli/Dm 

260  Jlratr=Rpr*Jl  lr-Rpi*Jl  li 

270  Jlrati=Rpr*Jlli+Rpi*Jllr 

280  OFF  ERROR 

290  CALL  Spherejnz(L- 1  ,U2r,U2i  J2rJ2i) 

300  CALL  Spherejnz(L,U2r,U2iJ21rJ21i) 

310  CALL  Spherehnz(L- 1  ,U2r,U2iT12rJi2i) 

320  CALL  Spherchnz(L,U2r,U2iJHL2  lr,H2 1  i) 

330  Nblr=L*J2r-(U2r*J21r-U2i*J21i) 

340  Nbli=L*J2i-(U2r*J21i+U2i*J21r) 

350  Nb2r=L-(U  lr*  Jlratr-U  li*Jlrati) 

360  Nb2i=-(Ulr*Jlrati+Uli*Jlratr) 

370  Nfacr=Mur*(Ep2r*J2r-Ep2i*J2i) 

380  Nfaci=Mur*(Ep2r*J2i+Ep2i*J2r) 

390  Topr=Nblr-(Nfacr*Nb2r-Nfaci*Nb2i) 
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400  Topi=Nbli-(Nfacr*Nb2i+Nfaci*Nb2r) 

410  Dblr=L*H2r-(U2r*H21r-U2i*H21i) 

420  Dbli=L*H2i-(U2r*H21i+U2i*H21r) 

430  Dfacr=Mur*(Ep2r*H2r-Ep2i*H2i) 

440  Dfaci=Mur*(Ep2r*H2i+Ep2i*H2r) 

450  Btmr=Db  lr-(Dfacr*Nb2r-Dfaci*Nb2i) 

460  Btmi=Dbli-(Dfacr*Nb2i+Dfaci*Nb2r) 

470  Bm=Btmr*Btmr+Btmi*Btmi 

480  Rbmr=Btmr/Bm 

490  Rbmi=-Btmi/Bm 

500  Ndma(L,  l)=-(Topr*Rbmr-Topi*Rbmi) 

510  Ndma(L,2)=-(Topr*Rbmi+Topi*Rbmr) 

520  NEXT  L 

530  SUBEND 

10  SUB  Ndeicon(Lmax,K2r,K2i,R?Ndea(*)) 

20  U2r=K2r*R 

30  U2i=K2i*R 

40  FOR  L=1  TO  Lmax 

50  CALL  Spherejnz(L- 1  ,U2r,U2i  J2r,J2i) 

60  CALL  Spherehnz(L- 1  ,U2r,U2i,H2r,H2i) 

70  Den=H2r*H2r+H2i*H2i 
80  Ih2r=H2r/Den 
90  Ih2i=-H2i/Den 
100  Ndea(L,  1  )=-(J2r*Ih2r-J2i*Ih2i) 

110  Ndea(L,2)=-(J2r*Ih2i+J2i*Ih2r) 

120  NEXT  L 
130  SUBEND 

140  SUB  Ndmicon(Lmax,K2r,K2i,R,Ndma(*)) 

150  U2r=K2r*R 

160  U2i=K2i*R 

170  FOR  L=1  TO  Lmax 

1 80  CALL  Spherejnz(L- 1  ,U2r,U2i  J2r,  J2i) 

190  CALL  Spherejnz(L,U2r,U2iJ21rJ21i) 

200  CALL  Spherehnz(L-LU2r,U2i>H2r>H2i) 

210  CALL  S  pherehnz(L,U  2r,U2i  di2 1  rJET2 1  i) 

220  A=U2r*H2 1  r-U  2i*H2 1  i 

230  B=U2r*H21i+U2i*H21r 

240  C=H2r*L-A 

250  D=H2i*L-B 

260  Den=C*C+D*D 

270  Ihr=C/Den 

280  Ihi=-D/Den 

290  A=U2r*J21r-U2i*J21i 

300  B=U2r*J21i+U2i*J21r 

310  C=J2r*L-A 

320  D=J2i*L-B 

330  Ndma(L,  1  )=-(C*Ihr-D*Ihi) 

340  Ndma(L,2)=-(C*Bii+D*Ihr) 

350  NEXT  L 
360  SUBEND 
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SERIES  REPRESENT ATIONS(v <  10, |z|  <  10) 

The  Bessel  functions  of  the  first  kind  7v(z)  can  be  constructed  from  the  series  representation 


Uz)  = 


{2)  kt0k\r(y  +  k  +  l). 


(Al) 


The  Hankel  function,  or  Bessel  function  of  the  third  kind  for  "outgoing  waves",  H^\z)  is 
related  to  /v(z)  and  the  Bessel  functions  of  the  second  kind,  Tv(z),  by 


W;i|(z)=J,(z)+iTv(2). 


(A2) 


For  non-integral  values  of  v,  Fv(z)  can  be  defined  by 

J y(z )  cos(vti)  — ■  J _y(Z  ) 
n(z)= — — 

but  for  integer  values  of  v,  it  is  necessary  to  use  the  more  complicated  expression 


(A3) 


(-0 


(A4) 


where  \y(/i)  is  given  by 


W  )®-7  .  ¥(«)  =  -JY+  i  T  (»S2),  (A5) 

k  =  lK 

and  yis  the  Euler  constant  0.577215664901532860606512. 

In  the  following  seven  subroutines,  the  functions /v(z)  and//v,}(z)  are  generated.  Use  of  these 
routines  has  been  restricted  to  cases  where  the  magnitudes  of  v  and  z  do  not  exceed  10.  For  larger 
values  of  v  and  z,  various  asymptotic  techniques  are  used.  The  discussion  of  these  begins  below, 
after  the  series  expansion  routines. 
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10  SUB  Jn(N,Zzr,ZziJnr,Jni,Err) 

20  !N  is  the  (Integer)  order. 

30  !Zzr  and  Zzi  are  die  real  and  imaginary  parts  of  the  argument. 

40  !  Jnr  and  Jni  are  the  real  and  imaginary  parts  of  the  Bessel  function. 

50  !Err  is  the  relative  truncation  error  in  the  truncation  of  the  series 

60  ! expansion  of  the  Bessel  function. 

70  M-0 

80  Zr=Zzr/2 

90  Zi=Zzi/2 

100  U2=Zi*Zi-Zr*Zr 

110  V2=-2*Zr*Zi 

120  IF  N=0  THEN 

130  Tr=l 

140  Ti=0 

150  Dr=Tr 

160  Di=Ti 

170  ELSE 

180  U=1 

190  V=0 

200  FOR  L=1  TO  N 

210  Ut=U*Zr-V*Zi 

220  Vt=U*Zi+V*Zr 

230  U=Ut/L 

240  V=Vt/L 

250  NEXT  L 

260  Tr=U 

270  Ti=V 

280  Dr=Tr 

290  Di=Ti 

300  END  IF 

310  M=M+1 

320  Ttr=(Tr*U2-Ti*V2)/M/(M+N) 

330  Tti=(Tr*V2+Ti*U2)/M/(M+N) 

340  Tr=Ttr 
350  Ti=Tti 
360  Dr=Dr+Tr 
370  Di=Di+Ti 

380  IF  (Tr*Tr+Ti*Ti)/(Dr*Dr+Di*Di)<Err*Err  THEN  400 

390  GOTO  310 

400  Jnr=Dr 

410  Jni=Di 

420  SUBEND 


10  SUB  Hln(N,Zzr^zi4ilnr,Hlni,EiT) 

20  !N  is  the  (Integer)  order 

30  !Zzr  and  Zzi  are  the  real  and  imaginary  parts  of  the  argument 
40  !Hlnr  and  Hlni  are  the  real  and  imaginary  parts  of  the  Bessel 
50  !  function.  Err  is  the  relative  truncation  error  in  the  truncation 

60  !of  the  series  expansion  of  the  Bessel  function. 

70  M=0 
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80  Zr=Zzr/2 
90  Zi=Zzi/2 
100  P1=0 

110  P2=0 

120  CALL  Prinlogz(Zr,Zi,LrJLi) 

130  Gam=.577215664901533 

140  Cr=PI-2*Li 
150  Cii=2*(Gam+Lr) 

160  U2=Zi*Zi-Zr*Zr 

170  V2=-2*Zr*Zi 

180  IF  N=0  THEN 

190  Sr=0 

200  Si=0 

210  Pwr=l 

220  Pwi=0 

230  Dr=Cr 

240  Di=Cii 

250  ELSE 

260  U=1 

270  V=0 

280  Nm=l 

290  FOR  L=1  TO  N 

300  Ut=U*Zr-V*Zi 

310  Vt=U*Zi+V*Zr 

320  U=Ut/L 

330  V=Vt/L 

340  P2=P2+1/L 

350  Nm=Nm*L 

360  NEXT  L 

370  Dn=U*U+V*V 

380  Uu=-V/Dn/Nm 

390  Vv=-U/Dn/Nm 

400  Nm=Nm/N 

410  Pwr=U 

420  Pwi=V 

430  Sr=Nm*Uu 

440  Si=Nm*Vv 

450  Dr=Cr*Pwr-(Cii-P2)*Pwi+Sr 

460  Di=Cr*Pwi+(Cii-P2)*Pwr+Si 

470  END  IF 

480  M=M+1 

490  P1=P1+1/M 

500  P2=P2+1/(N+M) 

510  Pwrt=(Pwr*U2-Pwi*V2)/M/(M+N) 

520  Pwit=(Pwr*V2+Pwi*U2)/M/(M+N) 

530  Pwr=Pwrt 

540  P^t=Pwit 

550  G=Cii-Pl-P2 

560  Tr=Pwr*Cr-Pwi*Ci 

570  Ti=Pwi*Cr+Pwr*Ci 

580  Dr=Dr+Tr 

590  Di=Di+Ti 
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600  IF  M<N  THEN 

610  Ssr=-(Sr*U2-Si*V2)/M/(N-M) 

620  Ssi=-(Sr*V2+Si*U2)/M/(N-M) 

630  Sr=Ssr 
640  Si=Ssi 
650  ELSE 
660  Sr=0 
670  Si=0 
680  END  IF 
690  Dr=Dr+Sr 
700  Di=Di+Si 

710  IF  (Tr*Tr+Ti*Ti)/(Dr*Dr+Di*Di)<Err*Err  THEN  730 

720  GOTO  480 

730  Hlnr=Dr/PI 

740  Hlni=Di/PI 

750  SUBEND 


10  SUB  Jnu(Nu,Zzr,Zzi  Jnur  JnuiJErr) 

20  !Nu  is  the  order 

30  !Zzr  and  Zzi  are  the  real  and  imaginary  parts  of  the  argument 
40  !  Jnur  and  Jnui  are  the  real  and  imaginary  parts  of  the  Bessel 

50  [function.  Err  is  the  relative  truncation  error  in  the  truncation 
60  !  of  the  series  expansion  of  the  Bessel  function. 

70  M=0 

80  Zr=Zzr/2 

90  Zi=Zzi/2 

100  U2=Zi*Zi-Zr*Zr 

110  V2=-2*Zr*Zi 

120  U=1 

130  V=0 

140  Dn-1 

150  Dr=l 

160  Di=0 

170  M=M+1 

180  Ut=U*U2-V*V2 

190  Vt=U*V2+V*U2 

200  U=Ut 

210  V=Vt 

220  Dn=Dn*(Nu+M)*M 
230  Tr=U/Dn 
240  Ti=V/Dn 
250  Dr=Dr+Tr 
260  Di=Di+Ti 

270  IF  ( ABS  (Tr)+ ABS  (Ti))/( ABS(Dr)+ AB  S  (Di))<Err  THEN  290 

280  GOTO  170 

290  Jnur=Dr 

300  Jnui=Di 

310  N=INT(Nu) 

320  Eps=Nu-N 

330  CALL  Garama(N  u^  1  ,G0) 
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340  CALL  Prinlogz(Zr,Zi,Lzr,Lzi) 
350  Mo=EXP(Eps*Lzr) 

360  Ph=Eps*Lzi 
370  Cr=Mo*COS(Ph) 

380  Ci=Mo*SIN(Ph) 

390  M=0 

400  Jr=(Cr*Jnur-Ci*Jnui)/G0 
410  Ji=(Cr*Jnui+Ci*Jnur)/G0 
420  U=1 
430  V=0 

440  IF  N<0  THEN  570 

450  IF  N=0  THEN  670 

460  U=Zr 

470  V=Zi 

480  M=M+1 

490  IF  M=N  THEN  670 

500  Ut=U*Zr-V*Zi 

510  Vt=U*Zi+V*Zr 

520  U=Ut 

530  V=Vt 

540  M=M+1 

550  IF  M=N  THEN  670 

560  GOTO  500 

570  Ut=U*Zr-V*Zi 

580  Vt=U*Zi+V*Zr 

590  U=Ut 

600  V=Vt 

610  M=M-1 

620  IF  M=N  THEN  640 

630  GOTO  570 

640  Dnm=U*U+V*V 

650  U=U/Dnm 

660  V=-V/Dnm 

670  Jnur=U*Jr-V*Ji 

680  Jnui=U*Ji+V*Jr 

690  SUBEND 


10  SUB  Hlnu(Nu,Zr,Zi,Hlnur,Hlnui^rr) 
20  ISeeJnu. 

30  CALL  Jnu(Nu,Zr^i  JrJiJErr) 

40  CALL  Jnu(-NuZr^i,JmrJmi£rr) 

50  C=COS(Nu*PI) 

60  S=SIN(Nu*PI) 

70  H 1  nur=Jr-(C*  Ji- Jmi)/S 
80  H 1  nui=Ji+(C*  Jr-Jmr)/S 
90  SUBEND 
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10  SUB  Prinlogz(Zr,Zi,LzrJLzi) 
20  Lzr=LOG(Zr*Zr+Zi*Zi)/2 
30  CALL  Atn2(Zr,ZifLzi) 

40  SUBEND 


10  SUB  Atn2(X,Y, Theta) 
20  IF  X=0  THEN 
30  IF  Y<0  THEN 
40  Theta=-PI/2 

50  ELSE 
60  Theta=PI/2 

70  END  IF 
80  GOTO  200 
90  END  IF 
100  Theta= ATN  (Y /X) 
110  IF  X>=0  THEN 
120  Theta=Theta 
130  ELSE 
140  IF  Y>=0  THEN 
150  Theta=Theta+PI 

160  ELSE 
170  Theta=Theta-PI 

180  END  IF 
190  END  IF 
200  SUBEND 


10  SUB  Gamma(X,G) 

20  OPTION  BASE  1 
30  ! 

40  ’.Computes  Gamma(x)  for  any  real  x( where  defined). 
50  IPrecision  is  IE-12 
60  ! 

70  DIM  A(20),Z(20) 

80  IF  X=1  THEN 
90  G=1 

100  GOTO  540 
110  END  IF 
120  IF  X>=1  THEN 
130  Xx=X 
140  Y=1 

150  Xx=Xx-l 

1 60  IF  Xx>=0  AND  Xx<=l  THEN  GOTO  280 

170  Y=Xx*Y 

180  Xx=Xx-l 

190  GOTO  160 

200  ELSE 

210  Xx=X-l 

220  Y=1 
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230  Xx=Xx+l 
240  Y=Y/Xx 

250  IF  Xx>=0  AND  Xx<=l  THEN  280 

260  GOTO  230 

270  END  IF 

280  A(l)=l 

290  A(2)=.57721566490 

300  A(3)=-.65587807 1 52 

310  A(4)=-.04200263503 

320  A(5)=.16653861 138 

330  A(6)=-.04219773456 

340  A(7)=-.00962197153 

350  A(8)=.00721 894325 

360  A(9)=-.001 16516759 

370  A(10)=-.0002 1524 167 

380  A(1 1)=.000 12805028 

390  A( 1 2)=-.000020 1 3485 

400  A(13)=-.00000125049 

410  A(14)=.000001 13303 

420  A(15)=-.00000020563 

430  A(16)=.00000000612 

440  A( 17)=.00000000500 

450  A(18)=-.OOOOOOOOU8 

460  A(19)=.00000000010 

470  A(20)= .00000000001 

480  Z(l)=Xx 

490  FOR  1=2  TO  20 

500  Z(I)=Xx*Z(I-l) 

510  NEXT  I 
520  Gi=DOT(A,Z) 

530  G=Y*Xx/Gi 
540  SUBEND 


ASYMPTOTIC  REPRESENTATIONS  (v  <  10,|z|  >  10) 


For  v  <  10  and  large  values  of  z,  the  following  asymptotic  expansions  are  useful 

Jv(z)  =  A/  — {P(y,  z)  cos  x  -  Q  (v,  z)  sin  x}  (A6) 

V  jiz 


and 


H?\z)  =  \[^-{P(v,z)+iQ(y,z)}eix  (A7) 

V  Jiz 

where  X  =  z“(iv+i)Jt  ^  11  =  4V2 ,  the  functions  P  (v,  z)  and  Q  (v,  z)  are  given  by 
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/>(V,z)- Xo(-l)  (2zf 

(A8) 

Q(V>Z)- tI(  1) 

(A9) 

and 

(v,0)=l  ,  (v,2)  =  (h-1)4^  ,  (v,4)  =  (|l-l)(tl-9)(n-25)fi^P  -  (A10) 
and 

(v,l)  =  H-l  ,  (v,3)  =  (n-l)0i-9)S^P  -  .  (All) 

The  next  four  subroutines  apply  these  expressions  to  the  calculation  of  the  Bessel  functions. 

10  SUB  Jnasy(N,Zr,Zi  Jnr  Jm,Trms) 

20  !Note:  Valid  for  arbitrary  real  N 

30  !Zr  and  Zi  real  and  imaginary  parts  of  argument 

40  !  Jnr  and  Jni  real  and  imaginary  parts  of  the  Bessel  function. 

50  ITrms  is  the  number  of  terms  retained  in  the  asymptotic  expansion. 

60  CALL  Pnz(N,Zr,Zi,Pnr,Pni,Trms) 

70  CALL  Qnz(N,Zr,Zi,Qnr,Qni,Trms) 

80  Chir=Zr-(2*N+l)*PI/4 
90  Chii=Zi 
100  E=EXP(Chii)/2 
110  Ei=  1/4/E 
120  C=COS(Chir) 

130  S=SIN(Chir) 

140  Cor=C*(E+Ei) 

150  Coi=-S*(E-Ei) 

160  Sir=S*(E+Ei) 

170  Sii=C*(E-Ei) 

180  Utlr=Pnr*Cor-Pni*Coi-Qnr*Sir+Qni*Sii 
190  UtIi=Pnr*Coi+Pni*Cor-Qnr*Sii-Qni*Sir 
200  Zm=SQR(Zr*Zr+Zi*Zi) 

210  Aa=SQR((Zm+Zr)/2) 

220  Bb=SQR((Zm-Zr)/2) 

230  Cc=Aa*Aa+Bb*Bb 

240  Dd=SQR(2/PI)/Cc 

250  Jnr=Dd*(Utlr*Aa+Utli*Bb) 

260  Jni=Dd*(-Utlr*Bb+Utli*Aa) 

270  SUBEND 


A-26 


NCSCTR  426-90 


10  SUB  Hlnasy(N,Zr,Zi,Hlnr,Hlni,Tnns) 

20  !Note:  Valid  for  arbitrary  real  N 

30  !Zr  and  Zi  real  and  imaginary  parts  of  argument 

40  !Hlnr  and  Hlni  real  and  imaginary  parts  of  the  Bessel  function. 

50  !Trms  is  the  number  of  terms  retained  in  the  asymptotic  expansion. 
60  CALL  Pnz(N^r,Zi,Pnr,Pni,Trms) 

70  CALL  Qnz(N,Zr£i,Qnr,Qni,Trms) 

80  Chir=Zr-(2*N+l)*PI/4 

90  Chii=Zi 

100  E=EXP(-Chii) 

110  C=COS(Chir) 

120  S=SIN(Chir) 

130  Utlr=(Pnr-Qni)*C-(Pni+Qnr)*S 
140  Utli=(Pnr-Qni)*S+(Pni+Qnr)*C 
150  Zm=SQR(Zr*Zr+Zi*Zi) 

160  Aa=SQR((Zm+Zr)/2) 

170  Bb=SGN(Zi)*SQR((Zm-Zr)/2) 

180  Cc=Aa*Aa+Bb*Bb 

190  Dd=E*SQR(2/PI)/Cc 

200  Hlnr=Dd*(Utlr*Aa+Utli*Bb) 

210  Hlni=Dd*(-Utlr*Bb+Utli*Aa) 

220  SUBEND 


10  SUB  Pnz(N,Zr,Zi,Pnr,Pni,Trms) 
20  !Used  with  Jnasy  and  Hlnasy. 

30  M=0 

40  Mu=4*N*N 
50  Tr=l 
60  Ti=0 
70  Dr=Tr 
80  Di=Ti 

90  U2=64*(Zr*Zr-Zi*Zi) 

100  V2=128*Zr*Zi 

110  Dn=U2*U2+V2*V2 

120  U=-U2/Dn 
130  V=V2/Dn 
140  M=M+2 

150  Ttr=Tr*U-TivV 
160  Tti=Tr*V+Ti*U 
170  Tr=Ttr/M/(M-l) 

180  Ti=Tti/M/(M-l) 

190  Fac=l 

200  FOR  1=1  TO  2*M-1  STEP  2 
210  Fac=Fac*(Mu-I*I) 

220  NEXT  I 

230  Dr=Dr+Tr*Fac 

240  Di=Di+Ti*Fac 

250  IF  M/2+1  <Trms  THEN  140 
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260  Pnr=Dr 
270  Pni=Di 
280  SUBEND 


10  SUB  Qnz(N,Zr,Zi,Qnr,Qni,Tniis) 
20  !Used  with  Jnasy  and  Hlnasy. 

30  M=1 

40  Mu=4*N*N 

50  Dnm=8*(Zr*Zr+Zi*Zi) 

60  Tr=Zr/Dnm 
70  Ti=-Zi/Dnm 
80  Dr=Tr*(Mu-l) 

90  Di=Ti*(Mu-l) 

100  U2=64*(Zr*Zr-Zi*Zi) 

110  V2=128*Zr*Zi 

120  Dn=U2*U2+V2*V2 

130  U=-U2/Dn 
140  V=V2/Dn 
150  M=M+2 

160  Ttr=Tr*U-Ti*V 
170  Tti=Tr*V+Ti*U 
180  Tr=Ttr/M/(M-l) 

190  Ti=Tti/M/(M-l) 

200  Fac=l 

210  FOR  1=1  TO  2*M-1  STEP  2 
220  Fac=Fac*(Mu-I*I) 

230  NEXT  I 

240  Dr=Dr+Tr*Fac 

250  Di=Di+Ti*Fac 

260  IF  (M+l)/2<Trms  THEN  150 

270  Qnr=Dr 

280  Qni=Di 

290  SUBEND 


UNIFORM  SYMPTOTIC  REPRESENTATION  (v >  10, | z|  >  10, v >\z\) 

Abramowitz  and  StegunA1  give  uniform  asymptotic  expansions  for  the  modified  Bessel 
functions  7v(z)  and  Ky(z).  These  functions  can  be  used  to  generate  the  functions  Jv(z)  and  H^Xz )  by 
means  of  the  identities 


and 


(A12) 


(A13) 
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where  the  argument  of  z  is  restricted  to  -n/2  <  arg  z  <  tc.  This  restriction  causes  no  problem,  since 
the  argument  of  z  in  electromagnetic  applications  is  limited  to  0  <,  arg  z  <  n/2. 

The  uniform  asymptotic  expansions  for  7v(z)  and  Ks(z)  are  given  by  (w  =  z/v) 


,,,  «vn  f, 

/v(z) — =====A  1  +  L  — j- 

V27tWl  +  wH  i=1  v* 


(A14) 


f,  ,  v  (~l)kUkit) 

'  h«+7l  *«  V* 


(A15) 


where 


1  1  w 
t  =  and  Tj  = — + In  ,  i  . 

r  1+t 


(A16) 


The  asymptotic  representation  of  the  product  of/v  and  77^  can  be  constructed  from  the  product 
of  their  asymptotic  representations,  and  this  will  help  avoid  overflow  and  underflow  problems  in 
computations.  The  product  is 


UMSXz’) =^e*,ni{zeTy{,'ef) 


(A17) 


where 


7v(z)/s:/z0~ 


,(vh-v'ti') 


2VwW(l+w2)(l  +  H''2)l  t=1  V* 


1+i^l|1+iH)Wr2 


*  =  l  v 


(A  18) 


The  functions  uk(t)  are  defined  by 


Uo(t)  =  1 


(A  19) 


«*+,(')  =  ^2d  -r2)u*(0+^  Jo  dx(l-^)uk( 


(A20) 


Abrcmcwitz  and  StegunA1  give  the  functions  uk  for  values  of£  =  0  thru  4  and  give  references 
for  £  =  5  and  6.  The  following  is  a  general  method  for  generating  uk  for  any  k.  Use  the  representation 
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uk(t)  —  ak0t  +ak'\t 


k+2*...+ak/+*  +  ...ak/+2k. 


Then  the  coefficients  are  given  by  the  following  recursion  formulas 


+  lpi 


a*+i.°  a*>c^2  +  8(Jk  +  l)] 

k+2i  1  1  I"  k  +  2(i-l)  _  J__  1 

2  +80fc  +  l+2/)  2  &(k  +  1+2/)  J 


and 


ak+i,t+ 1 


k+2k  5 

2  +  8(£  +  l+2{£  +  l}) 


where  i  =  1,2, 


10  SUB  Jfiuniasym(Nu,Zr,Zi,Jr,Ji,Trins) 

20  !Nu  is  the  order. 

30  !Zr  and  Zi  real  and  imaginary  parts  of  argument 

40  !  Jr  and  Ji  real  and  imaginary  parts  of  the  Bessel  function. 

50  !Trms  is  the  number  of  terms  retained  in  the  asymptotic  expansion. 
60  DIM  A(10,10),U(10,2),T(30,2) 

70  REDIM  A(Trms,Trms),U(Trms,2),T(3*Trms,2) 

80  Wr=Zi/Nu 

90  Wi=-Zr/Nu 

100  Ur=l+Wr*Wr-Wi*Wi 

110  Ui=2*Wr*Wi 

120  Um=SQR(Ur*Ur+Ui*Ui) 

130  Radr=SQR((Um+Ur)/2) 

140  Radi=SGN(Ui)*SQR((Um-Ur)/2) 

150  Radm=Radi*Radi+Radr*Radr 
160  Sradm=SQR(Radm) 

170  Radrtr=SQR((Sradm+Radr)/2) 

180  Radrti=SGN(Radi)*SQR((Sradm-Radr)/2) 

190  Radrtm=Radrtr  *  Radrtr +Radrti  *  Radrti 

200  Iraditr=Radrtr/Radrtm 

210  Iradrti=-Radrti/Radrtm 

220  Tr=Radr/Radm 

230  Ti=-Radi/Radm 

240  Dnr=l+Radr 

250  Dni=Radi 

260  Dnm=Dnr*Dnr+Dni*Dni 

270  Nmr=Dnr/Dnm 

280  Nmi=-Dni/Dnm 


(A21) 

(A22) 

(A23) 

(A24) 
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250  Dni=Radi 

260  Dnm=Dnr*Dnr+Dni*Dni 

270  Nmr=Dnr/Dnm 

280  Nmi=-Dni/Dnm 

290  Argr=Wr*Nmr-Wi*Nmi 

300  Argi=Wr*Nmi+Wi*Nmr 

310  CALL  Prinlogz(Argr,  ArgiJ.gr,  Lgi) 

320  Etar=Radr+Lgr 

330  Etai=Radi+Lgi 

340  CALL  Ukcoefs(Trms,A(*)) 

350  CALL  Uk(Tr,Ti,Trms  A(*),T(*),U(*)) 

360  Dumr=0 

370  Dumi=0 

380  Nufac=Nu 

390  FOR  1=1  TO  Trms 

400  Nufac=Nufac/Nu 

410  Dumr=Dumr+U(I,l)*Nufac 

420  Dumi=Dumi+U(I,2)*Nufac 

430  NEXT  I 

440  E=EXP(Nu*Etar) 

450  Co=COS  (N  u*Etai) 

460  Si=SIN(Nu*Etai) 

470  Fac=l/SQR(2*PI*Nu) 

480  Facr=Fac*E*(Co*Iradrtr-Si*Iradrti) 
490  Faci=Fac*E*(Co*Iradm+Si*Iradrtr) 
500  Inur==Facr*Dumr-Faci*Dumi 
510  Inui=Facr*Dumi+Faci*Dumr 
520  Cn=COS  (N  u*PI/2) 

530  Sn=SIN(Nu;tPI/2) 

540  Jr=Cn*Inur-Sn*Inui 
550  Ji=Cn*Inui+Sn*Inur 
560  SUBEND 


10  SUB  Hfkuniasyin(Nu,Zr  ,Zi,H  lr,H  li,Trms) 

20  !Nu  is  the  order. 

30  !Zr  and  Zi  real  and  imaginary  parts  of  argument 
40  !Hlr  and  Hli  real  and  imaginaj7  parts  of  the  Bessel  function. 

50  !Trms  is  the  number  of  terms  retained  in  the  asymptotic  expansion. 

60  DIM  A(10,10),U(10,2),T(30,2) 

70  REDIM  A(Trms,Trms),U(Trms^),T(3*Trms,2) 

80  Wr=Zi/Nu 

90  Wi=-Zr/Nu 

100  Ur=l+Wr*Wr-Wi*Wi 

110  Ui=2*Wr*Wi 

120  Um=SQR(Ur*Ur+Ui*Ui) 

130  Radr=SQR((Um+Ur)/2) 

140  Radi=SGN(Ui)*SQR((Um-Ur)/2) 

150  Radm=Radi  *  Radi+Radr  *Radr 
160  Sradm=SQR(Radm) 

170  Radrtr=SQR((Sradm+Radr)/2) 
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180  Radrti=SGN(Radi)*SQR((Sradm-Radr)/2) 

190  Raditm=Raditr*Radrtr+Radrti*Radrti 

200  Iradrtr=Radrtr/Radrtm 

210  Iradrti=-Radrti/Radrtm 

220  Tr=Radr/Radm 

230  Ti=-Radi/Radm 

240  Dnr=l+Radr 

250  Dni=Radi 

260  Dnm=Dnr*Dnr+Dni*Dni 

270  Nmr=Dnr/Dnm 

280  Nmi=-Dni/Dnm 

290  Argr=Wr*Nmr-Wi*Nmi 

300  Argi=Wr*Nmi+Wi*Nmr 

310  CALL  Prinlogz(Argr,Argi,Lgr,Lgi) 

320  Etar=Radr+Lgr 
330  Etai=Radi+Lgi 
340  CALL  Ukcoefs(Trms,A(*)) 

350  CALL  Uk(Tr,Ti»TrmsA(*).T(*),U(*)) 

360  Dumr=0 

370  Dumi=0 

380  Nufac=-Nu 

390  FOR  1=1  TO  Tims 

400  Nufac=-Nufac/Nu 

410  Dumr=Dumr+U(I,l)*Nufac 

420  Dumi=Dumi+U(I,2)*Nufac 

430  NEXT  I 

440  E=EXP(-N  u*Etar) 

450  Co=COS(Nu*Etai) 

460  Si=-SIN(Nu*Etai) 

470  Fac=SQR(PI/2/Nu) 

480  Facr=Fac*E*(Co*Iradrtr-Si*Iradrti) 

490  Faci=Fac*E*(Co*Iradrti+Si*Iradrtr) 

500  Knur=Facr*Dumr-Faci*Dumi 
510  Knui=Facr*Dumi+Faci*Dumr 
520  Cn=COS(Nu*PI/2) 

530  Sn=-SIN(Nu*PI/2) 

540  Jr=Cn*Knur-Sn*Knui 
550  Ji=Cn*Knui+Sn*Knur 
560  Hlr=2*Ji/PI 
570  Hli=-2*Jr/PI 
580  SUBEND 


10  SUB  Jhlproduni(Nu,Niip,Zr,Zi,Zpr,Zpi, Prodr, Prodi,Tnns) 

20  !Nu  and  Nup  are  the  orders  of  the  two  Bessel  functions  in  the 
product 

30  !Zr  and  Zi  are  the  real  and  imaginary  parts  of  the  argument  of  one 
40  !  factor  and  Zpr  and  Zpi  are  those  for  the  other  factor. 

50  !  Prodr  and  Prodi  are  die  real  and  imaginary  parts  of  the  Bessel 

60  !  function  product 

70  ITrms  is  the  number  of  terms  retained  in  the  asymptotic  expansion. 
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80  DIM  A(10,10),U(10,2),T(30,2) 

90  REDIM  A(Trms,Tnns),U(Trms,2),T(3*Trms,2) 

100  Wr=Zi/Nu 

110  Wi=-Zr/Nu 

120  Ur=l+Wr*Wr-Wi*Wi 

130  Ui=2*Wr*Wi 

140  Um=SQR(Ur*Ur+Ui*Ui) 

150  Radr=SQR((Um+Ur)/2) 

160  Radi=SGN(Ui)*SQR((Um-Ur)/2) 

170  Radm=Radi  *  Radi+Radr*  Radr 
180  Sradm=SQR(Radm) 

190  Radrtr=SQR((Sradm+Radr)/2) 

200  Radrti=SGN(Radi)*SQR((Sradm-Radr)/2) 

210  Radrtm=Radrtr  *  Radrtr+Radrti  *  Radrti 

220  Iradrtr=Radrtr/Radrtm 

230  Iradrti=-  Radrti/Radrtm 

240  Tr=Radr/Radm 

250  Ti=-Radi/Radm 

260  Dnr=l+Radr 

270  Dni=Radi 

280  Dnm=Dnr*Dnr+Dni*Dni 

290  Nmr=Dnr/Dnm 

300  Nmi=-Dni/Dnm 

310  Argr=Wr*Nmr-Wi*Nmi 

320  Argi=Wr*Nmi+Wi*Nmr 

330  CALL  Prinlogz(Argr,ArgiLgrJLgi) 

340  Etar=Radr+Lgr 

350  Etai=Radi+Lgi 

360  CALL  Ukcoefs(Trms^(*)) 

370  CALL  Uk(Tr,Ti,Trms,A(*),T(*),U(*)) 

380  Dumr=0 

390  Dumi=0 

400  Nufac=Nu 

410  FOR  1=1  TO  Trms 

420  Nufac=Nufac/Nu 

430  Dumr=Dumr+U(I,l)*Nufac 

440  Dumi=Dumi+U(I,2)*Nufac 

450  NEXT  I 

460  MAT  U=  (0) 

470  Wr=Zpi/Nup 

480  Wi=-Zpr/Nup 

490  Ur=l+Wr*Wr-Wi*Wi 

500  Ui=2*Wr*Wi 

510  Um=SQR(Ur*Ur+Ui*Ui) 

520  Radr=SQR((Um+Ur)/2) 

530  Radi=SGN(Ui)*SQR((Um-Ur)/2) 

540  Radm=Radi*Radi+Radr*Radr 
550  Sradm=SQR(Radm) 

560  Radrtr=SQR((Sradm+Radr)/2) 

570  Radrti=SGN  (Radi)*  SQR((Sradm-Radr)/2) 
580  Radrtm=Radxtr*Radrtr+Radrti*Radrti 
590  Iradrtpr=Radrtr/Radrtm 
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600  Iradrtpi=-Radrti/Radrtm 

610  Tr=Radr/Radm 

620  Ti=-Radi/Radm 

630  Dnr=l+Radr 

640  Dni=Radi 

650  Dnm=Dnr*Dnr+Dni*Dni 

660  Nmr=Dnr/Dnm 

670  Nmi=-Dni/Dnm 

680  Argr=Wr*Nmr-Wi*Nmi 

690  Argi=Wr*Nmi+Wi*Nmr 

700  CALL  Prinlogz(Argr,Argi,Lgr,Lgi) 

710  Etapr=Radr+Lgr 
720  Etapi=Radi+Lgi 
730  CALL  Uk(Tr,Ti,Trms,A(*),T(*),U(*)) 

740  Dumpr=0 

750  Dumpi=0 

760  Nufac=-Nup 

770  FOR  1=1  TO  Trms 

780  Nufac=-Nufac/Nup 

790  Dumpr=Dumpr+U(I,l)*Nufac 

800  Dumpi=Dumpi+U(I,2)*iNufac 

810  NEXTI 

820  Sumr=Dumr*Dumpr-Dumi*Dumpi 

830  Sumi=Dumr*Dumpi+Dumi*Dumpr 

840  Radsr=Iradrtr*Iradrtpr-Iradrti*Iradrtpi 

850  Radsi=Iradrtr*Iradrtpi+Iradrti*Iradrtpr 

860  Dumprodr=(Sumr*Radsr-Sumi*Radsi)/(2*SQR(Nu*Nup)) 

870  Dumprodi=(Sumr*Radsi+Sumi*Radsr)/(2*SQR(Nu*Nup)) 

880  Etargr=Nu*Etar-Nup*Etapr 

890  Etargi=Nu*Etai-Nup*Etapi 

900  E=EXP(Etargr) 

910  Co=COS(Etargi) 

920  Si=SIN(Etargi) 

930  Ikprodr=E*(Dumprodr*Co-Dumprodi*Si) 

940  Ikprodi=E*(Dumprodr*Si+Dumpxodi*Co) 

950  Cfr=2*SIN((Nu-Nup)*PI/2)/PI 
960  Cfi=-2*COS((Nu-Nup)*PI/2)/PI 
970  I>rodr=Ikprcxlr*Cfr-Ikpr(xii*Cfi 
980  Prodi=Ikprodr*Cfi+Dqprodi*Cfr 
990  SUBEND 


10  SUB  Uk(Tr,Ti,K,A(*),T(*),Uk(*)) 

20  IFunctions  used  in  the  uniform  asymptotic  expansion  of  the  Bessel 
30  '.functions  for  large  orders.  A(*)  is  a  matrix  generated  by  the 
40  !  program  Ukcoefs.  T(*)  and  Uk(*)  are  output  to  the  uniform  asymp- 

50  Itotic  Bessel  function  routines. 

60  T(l,l)=Tr 

70  T(l,2)=Ti 

80  FOR  1=2  TO  3*K 

90  T(1, 1  )=T(I- 1 , 1  )*Tr-T(I- 1 ,2)*Ti 
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100  T(I,2)=T(I- 1 ,  l)*Ti+T(I- 1 ,2)*Tr 

110  NEXT  I 

120  Uk(l,l)=l 

130  Uk(l,2)=0 

140  FOR  1=2  TO  K 

150  FOR  J=1  TO  I 

160  Uk(1, 1  )=Uk(1, 1  )+A(IrT)*T(2*J+I-3, 1 ) 

170  Uk(I,2)=Uk(I,2)+A(I,J)*T(2*J+I-3,2) 

180  NEXT  J 
190  NEXT  I 
200  SUBEND 


10  SUB  Ukcoefs(K,A(*)) 

20  !K  is  the  maximum  index  value.  A(*)  is  output  to  Uk. 

30  A(l,l)=l 

40  FOR  L=0  TO  K-2 

50  A(L+2, 1  )= A(L+ 1 , 1  )*  (L/2+ l/8/(L+ 1 )) 

60  A(L+2,L+2)=- A(L+ 1  ,L+ 1 )  *  (3  *L/2+5/24/  (L+ 1 )) 

70  IF  L>=1  THEN 
80  FORM=l  TOL 

90  A(L+2,M+l)=A(L+l,M+l)*((L+2*M)/2+l/(8*(L+l+2*M)» 

100  A(L+2,M+l)=A(L+2,M+l)-A(L+l>l)!fc((L+2*(M-l))/2+5/(8*(L+l+2*M))) 

110  NEXT  M 

120  END  IF 

130  NEXT  L 

140  SUBEND 
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